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CHAPTER 1 


INTRODUCTION 


Complex structural systems are most often modeled for analysis as 
assemblies of discrete structural components. The most generally ap- 
plicable discretization approach is the finite element method (FEM). In 
this method It Is often necessary to divide the structural model Into a 
very large number of elements In order to accurately evaluate displace- 
ments, strains, and stresses. As the number of elements Increases, the 
number of degrees of freedom (DOF) in the model can easily exceed the 
capacity of many present-day computer facilities (both hardware and 
software) or can make the solution of the large order matrix equations 
prohibitively expensive [ 1 . 13 . This problem becomes particularly acute 
In nonlinear analysis. The Iterative nature of nonlinear analysis re- 
quires that the matrix equations be solved repetitively, thus com- 
pounding the computational expense. 

The application of time-dependent (dynamic) loading to complex 
structural models imposes additional difficulty on current structural 
software systems. A linear dynamic analysis can be orders of magnitude 
more expensive than a static analysis of the same model. When nonlinear 
response Is also considered, the computational effort required for 
analysis can quickly make the solution infeasible. 
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Recent advances in computer arch i tecture, primarily in brute com- 
putational speed, have somewhat relieved the problems of excessively 

large and costly analyses. It is not expected, however, that new 
developments in hardware will keep pace with the growing demands, in 
terms of model size and complexity, made by structural analysts. Nor 
will the very expensive super-computers become widely available. 
Moreover, it appears unlikely that orders of magnitude increases in 
high-speed memory capacity and data transmission rates (which analysts 
have come to expect every few years) can be extracted from current tech- 
nology. If any short term relief Is to come, it must derive from im- 
provements in the structural analysis software. Such improvements lie 
In more efficient use of existing hardware and in improved structural 
modeling techniques. The focus of this work is on the improvement of 
structural modeling techniques. 

One procedure that is used successfully in static analysis is mul- 
tilevel substructured modeling Cl .23. This approach allows the various 
major structural units, or substructures, to be treated independently 
prior to final assembly. With the use of condensation techniques, a 
reduction in total model size can be achieved while exactly retaining 
the original system characteristics for static analysis. Substructuring 
techniques also find broad applicability to the various types of com- 
puter hardware used by engineers today. Efficient use of both main- 
frames and virtual memory minicomputers with either serial or pipeline 
processors has been demonstrated Ci .33. In view of independent sub- 
structures, the possible adaptation of the software to a system of in- 
dependently operating processors (or computers) under the logical con- 
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fro I ot a single machine becomes quite attractive. 

Many investigators have presented extensions of the substructur i ng 
approach from static analysis to dynamic analysis of finite element 
models. However, these efforts are limited to one level of substruc- 
tures. No attempt has yet been made to formulate and implement mu I- 
' 1 1 1 eve I substructuring for dynamic analysis In a general purpose FEM 
system. Thus, the techniques are not proven effective from the prac- 
tical viewpoint of large scale structural analysis. The need exists for 
a comprehensive dynamic analysis system capable of processing multilevel 
substructured models and Incorporating nonlinear response Into the solu- 
tion. 


1 .2 Obiecti ves aM Scope 


The objectives of this work are fourfold: 

1. To review the state-of-the-art of the multilevel substructure 
methodology and modeling procedures for the static analysis of 
complex structures by the FEM. Included is a presentation of 
the design and implementation of the required software and an 
illustration of the modeling technique by way of example 
problems. 

2. To review the analytical formulations and computational 
procedures available for the analysis of complex structural 
systems subjected to time- dependent loads and capable of 
linear or nonlinear response. Emphasis is placed on methods 
for reducing the size of the finite element model for dynamic 
analysis. Also studied are efgenproblem solution procedures, 
solution of the equations of motion, and formulations for 
tracking nonlinear response. 

3. To identify the reduction and computational procedures most 
suitable for incorporation into a general FEM software system 
with multilevel substructur ing capabilities. 

4. To discuss the Implementation considerations of the above 
relative to a multilevel substructured environment. 
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The discussion of dynamic analysis methods is based on a review and 
i nterpretation of the open I iterature. While the authors, along with 
several other researchers, have experience in large scale linear and 
nonlinear static analysis with substructur ing, large scale dynamic 
analysis has been performed using only the simplest of techniques. 

The remainder of this report is divided Into chapters which iden- 
tify the major topics covered. Chapter 2 is a presentation of the 
methodology and implementation procedures of multilevel substructur ing 
In a general purpose software system. Methods for reducing the order of 
the coefficient matrices in the differential equations of motion are 
reviewed and evaluated in Chapter 3. Chapter 4 considers the various 
computational algorithms required In dynamic analysis. Included are 
eigenproblem solution methods, procedures for solving the differential 
equations of motion, and selected minor topics. The nonlinear continuum 
mechanics equations for finite deformation, cast in matrix form, are 
presented In Chapter 5. Exact forms for both the Total Lagrangian and 
the Updated Lagrangian approaches are described in addition to the 
finite deformation theories of elasto-plasticity. Chapter 6 presents 
specific forms of the finite element matrices for the Total and Updated 
Lagrangian approaches. Details of the transient solution procedure 
based on an implicit integration operator with the effects of substruc- 
turing are discussed. Chapter 7 describes the input language designed 
to provide a convenient user Interface with the application software. A 
summary of work performed thus far Is presented In Chapter 8 along with 
proposed future activities necessary for successful dynamic analysis of 
multilevel substructured models. 
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1 .3 Notat’i on 

Most of the notation used In the following chapters Is defined as 
It is Introduced. However, the following list may prove to be a useful 
reference. Consistency is maintained whenever possible and exceptions 
are noted In the text. The complexity of the notation used in Chapters 
5 and 6 warrants special discussion in those chapters. 


{ } 

[ ] 

pQ) 


{q} 

{u}, {u}, { u} 


{U m }, {u 5 } 


Cx j3 

K.E. 


S.E. 

DO 

CB], [G] 


Cc(t)] 

[d: 

DO, CM] 


a vector 
a matrix 

characteristic polynomial 
global coordinate system 
local coordinate system 

vector of generalized displacement coord i nates 
displacement, velocity, and acceleration vectors In 
geometric coordinates 

displacement vectors for master and slave DOF 
matrix of trial mode shapes for Iteration " j " 
kinetic energy for a structure 
strain energy for a structure 
a tr I diagonal matrix 

an Interaction matrix used in simultaneous Iteration 
methods 

time dependent damping matrix 
substructure dynamic stiffness matrix 
stiffness and mass matrices or submatrices 
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CKq3> CMq 3 Guyan reduced stiffness and mass matrices 
CKp], [Mp] f i xed- i nterf ace substructure reduced stiffness and mass 

matr i ces 

[fcj k ], [M^] submatrices of generalized stiffness and mass for sub- 
structure "i" 

* * 

[K], [M] stiffness and mass matrices assembled from lower level 

substructures 

x mm 

[K ], [M ] stiffness and mass submatrices for the system master DOF 

within [Kl] and CM] 


CL] 


CM(co)] 

CP], CQ] 

(PCt)} 

CR] 

Cs] 

Ct g ], C4> c ] 
CTJ 

w 2 

C u 2 ] 


x, DJ 

M 

C$ n J 


lower triangular matrix of Choleski factors of a stiffness 
matrix 

dynamic mass matrix 
orthogonal transformation matrices 
time dependent load vector 
an upper triangular matrix 

symmetric coefficient matrix for an eigenproblem in stan- 
dard form 

Guyan transformation matrix (static constraint modes) 

frequency dependent transformation matrix 

square of the substructure natural frequency for mode "l w 

diagonal matrix of substructure natural frequencies 

squared 

unknown system frequency 

an eigenvalue, diagonal matrix of eigenvalues 
matrix of substructure normal modes of vibration 
set of retained substructure normal modes 
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CHAPTER 2 


MULTILEVEL SUBSTRUCTURE STATIC ANALYSIS 


Complex structures frequently consist of repeated. Identical compo- 
nents. This may be dictated by economic, construction, or symmetry con- 
straints. The boundaries between components (either real or artificial) 
partition a complex structure Into a natural system of simpler substruc- 
tures. Each substructure may In turn be partitioned to exploit ad- 
ditional repetition. The associated flntte element model generation 
process Is considerably simplified through the repeated use of 
previously defined components. In many cases, the computational costs 
of analysis are reduced accordingly. This concept of substructured 
modeling has been termed the H su per element” technique In the literature 
due to the similarity of the substructure merging process with that used 
to connect finite elements. The term "user-defined" has also been em- 
ployed to distinguish analyst specified substructuring from automatic 
partitioning of the equilibrium equations. 

The structural frame of an aircraft provides the classic example to 
Illustrate the concepts and advantages of a substructured analysis. In- 
dependent design groups develop the Individual substructures, for ex- 
ample: the wing assembly, fusealage sections, and vertical stabilizer. 
The substructures interface through relatively small boundaries (In 
terms of the number of nodes). Even with such "f Irst- level " substruc- 
turing, the number of nodes and elements may be too large for efficient 
processing of the substructures. The same substructuring process. In 
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theory, can be repeated within each first- 1 eve I substructure to In- 
troduce second, third, fourth, ... level substructures. The engine 
nacelles, shrouds, flaps, ribs, and skin panel sections within a wing 
assembly may comprise substructures five or six levels removed from the 
"highest” level structure that represents the complete airframe. The 
•conceptual organization for this type of structural mode! parallels that 
of a tree. The tree has a single root node (the highest level struc- 
ture). Any number of substructure levels may be defined below the root 
node. No theoretical I tmit exists on the number of branches that enter 
a node (substructure) at level "l” from level "i-1". All terminal 
nodes of the tree are individual finite elements. 

Substructure techniques have been utilized extensively by the Nor- 
weigan ship building Industry £2.33 to construct and analyse finite ele- 
ment models of oil tankers. Repeated bulkheads and common stiffener ar- 
rangements in ships are well suited for substructur I ng. Without mul- 
tilevel substructures, typical analyses would Involve 100-150 thousand 
degrees of freedom. Problems of this size remain Impractical to solve 
despite the availability of super-computers. Both the airframe and ship 
building examples clearly demonstrate the usefulness of multi-level sub- 
structuring to support practical analyses. 

In the context of linear, static analysis, a substructured model 
yields the same solution (to within round-off errors) as a "standard" 
model that considers the structure as a single collection of nodes and 
elements. The solution equivalence remains valid for static, nonlinear 
analysis provided the substructures experiencing nonlinear behavior are 
included in the iterative solution. Normally, regions of a structure 



that remain linear are substructured and eliminated via static condensa- 
tion from the Iterative solution. For example, it Is a simple matter to 
anticipate the plastic deformation that occurs near a stress concentra- 
tion. Portions of the structure removed from the stress concentration 
are defined as linear substructures and condensed. Nonlinear finite 
•elements and reduced linear substructures comprise the highest level 
structure for the iterative solution. Linear substructures simply 
provide elastic restraint on the nonlinear region. In such cases, the 
standard and substructured models again yield Identical solutions. 
However, the substructured model generally requires much less com- 
putational' effort as a consequences of its reduced size. 

This chapter provides detailed background information on substruc- 
tured modeling and solution procedures for static analysis. The advan- 
tages of substructured analysis relative to the standard procedure are 
first described. The literature concerned with static substructuring Is 
reviewed. Various techniques adopted to address the user-software in- 
terface are described and a simple example is presented to Illustrate 
the techniques used In the POLO-FINITE system. Computational and 
software Issues that arise in programs for general substructured 
analysis are surveyed. The chapter concludes with a study of two ex- 
ample problems that illustrate the computational savings possible with 
substructured models for linear and nonlinear static analysis. 
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2 .2 Substructured vs Standard Mode 1 s 


Compared with a standard modeling and solution procedure, a multi- 
level substructured approach offers a number of advantages. These i n- 
cl ude: 


1. Input data requirements are reduced. Geometric and topo logic 
descriptions of a substructure are specified only once. When 
the same substructure appears repeatedly at higher levels, in- 
put data that must be provided by the user is significantly 
reduced. 

2. Tne impact of design changes on reanalysis costs can be 
reduced. Stiffness matrices and loading vectors for only the 
modified substructures are computed during reanalysis. 

3. Models may be generated independently. Because substructures 
have clearly defined interfaces, the design and modeling 
groups may work almost independently. Element numbering, node 
numbering, and load case naming schemes are usually Indepen- 
dent across substructures which simplifies model generation. 

4. Isolated substructures may be Independently verified. Each 
substructure may be constrained, loaded, and analyzed to study 
the response of the isolated model. Checkout analyses are 
usually inexpensive relative to the cost of a complete struc- 
ture analysis. The higher execution priority assigned to jobs 
that request fewer machine resources also decreases the total 
analysis time. A complete analysis performed in smaller seg- 
ments frequently costs less, and requires less residency time, 
than a comparable standard analysis performed In a single ex- 
ecution. 

5. Identical substructures may be used repeatedly. This Is the 
most often cited advantage of a substructured model. Sub- 
structure quantities, for example the reduced stiffness 
matrix, are computed but once and used repeatedly to form 
other structure stiffness matrices. The degree of com- 
putational savings Is highly problem dependent with cost 
reductions in the range of 2 to 100 having been reported 
C2 .2]. The level of savings in nonlinear analysis depends on 
the degree of size reduction and the frequency of tangent 
stiffness updates. 

6. Numerical conditioning problems are often reduced. Certain 
types of numerical problems may be remedied through the use of 
a substructured modeling and solution procedure. A common 
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situation involves the joining of a very stiff substructure to 
a comparatively flexible one which results in a loss of preci- 
sion during stiffness assembly (a very large number overwhelms 
the much smaller term). Condensation of a very stiff sub- 
structure frequently reduces the magnitude of stiffness coef- 
ficients for the remaining nodes to a level comparable with 
those of the adjacent, more flexible component. The loss of 
precision in the important stiffness assembly process can thus 
be minimized without extended precision arithmetic. 

7. Exposure to machine failure is reduced. The solution ot a 
structure with a very large number of DOF may require long 
residency times on multipurpose computer systems. During 
this time, the executing program is exposed to the possibility 
of a machine failure that would require restarting the 
analysis. Procedures have been developed, termed checkpoint 
restart, that save snapshots of the program status on disk 
files at specified intervals. This process often requires ex- 
tensive machine dependent coding which restricts the software 
portability. Substructured models provide a more natural 
solution to the failure protection problem. Substructures are 
processed In a logically independent, sequential manner during 
execution. Natural breakpoints occur between each substruc- 
ture at which the the execution may be terminated and the 
databases saved on tape. If a failure occurs while the next 
substructure is being processed, the databases on tape are 
simply restored to disk and the analysis resumed. 


These advantages of a substructured approach are equally applicable in 
all types of analyses — static, dynamic, linear, and nonlinear. 


2.3 


Sur v . e .y 


Although the concepts and equations of substructured analysis are 
generally well known for static analysis, relatively few papers on the 
subject have appeared In the literature. A contributing factor to this 
apparent lack of interest is the enormous software complexity required 
to support substructuring, coupled with the limited availability ot 
software to researchers. For static linear and nonlinear analysis, the 
governing equations are straightforward and simple to derive. The few 
researchers who have examined substructuring have focused on improving 
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the computational efficiency and on the computer implementation 
problems. Because only approximate substructured solutions are feasible 
in dynamic analysis, the literature concerned with improving the al- 
gorithms continues to expand and is reviewed in the next chapter. This 
section reviews previous studies that addressed static, linear and non- 
linear substructured analyses. 

In the early 1960s, Przemteniecki [2.15] presented the first com- 
prehensive formulation for substructured analysis using the conventional 
displacement method. Taig [2.17] described an attempt to implement 
these procedures for general analysis. These two early studies viewed 
substructuring as automatic partitioning of the equilibrium equations to 
overcome computer size limitations. Interest In substructured analysis 
declined during the late sixties when commercial programs using sparse 
matrix techniques became operational on third generation computers. The 
non-computational aspects of analysis. Including model generation and 
technical coordination, replaced computer limitations as the major 
problem areas. Substructure techniques regained appeal as an approach 
to solve these problems but in the form of analyst defined, rather than 
automatic partitioning. Thus the term "user-defined" substructuring was 
coined. In this same period, Williams [2.18] showed through operation 
counts that equation solving with sparse matrix techniques can never be 
more computationally efficient than a substructured solution, when the 
substructure arrangement is suitably defined. 
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Furuike [2.53 described a software system capable of processing 
three levels of substructures. The root and the next two levels of the 
tree may contain only substructures. Fourth level structures consist 
entirely of finite elements. The analyst supplies, through input data, 
the order of substructure processing and the proper sequencing of nodal 
DOF across substructure boundaries. Substructure stiffness condensation 
is performed with the inefficient Inversion algorithm described in Sec- 
tion 2.5.1. The author presents several example solutions with and 
without substructuring that demonstrate the numerical equivalence of the 
resuits. Unfortunately, no comparisons of computer resources are 
provided. 

Egelan and Araldsen [2.33 briefly surveyed the substructur i ng 
capabilities of the SESAM-69 program. No details of the user- interface 
for describing the substructure connectivity and orientation are 
provided. SESAM supports multilevel substructured models. Finite ele- 
ments and substructures cannot be mixed at the same level in SESAM. 
This restriction Is imposed by the adoption of separate assembly proces- 
sors for structures composed of only elements and those composed of only 
substructures. The classical Inversion algorithm for stiffness conden- 
sation Is Implied In the survey. 

Descriptions of substructuring capabilities In NASTRAN are given by 
MacNeal and McCormick [2.113 and in ASKA by Schrem [2.163. Each system 
Initially supported only modest substructuring; the major emphasis being 
directed toward efficient sparse matrix techniques. Proprietary ver- 
sions of NASTRAN have been expanded to support more comprehensive sub- 
structuring, Including some techniques for dynamic analysis beyond Guyan 
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reduction. Software details are not publicly available. To perform a 
substructured analysis, users must write "driver" programs In a special 
language to manipulate disk files and to control the execution of 
processing modules. Completely automated solution procedures for sub- 
structured analysis are not available. 

More recently, Peterson and Popov C2.143 addressed the com- 
putational penalty that occurs with stiffness matrix rearrangement prior 
to condensation. They propose a scheme to el Iminate nodal freedoms at 
the element level before assembly of the condensed stiffness. The tech- 
nique is promoted as more efficient than the conventional method of 
rearrangement although no comparisons of computer time are stated for 
the example problems. 

Lopez [2.93 presented POLO-FINITE as the first major system to sup- 
port multi-level substructuring as a natural approach to finite element 
modeling. Major advances incorporated In POLO-FINITE include a very 
simple substructure definition process and fully automated solution 
procedures. The dependencies between substructures in the hierarchy are 
determined completely by the system from basic Input. Any number of 
related and/or independent hierarchies may be defined within a single 
database, with up to 20 levels of substructuring permitted in each 
hierarchy. Substructures and finite elements may be mixed at any level; 
the system processors treat finite elements and substructures Iden- 
tically. The extensive system logic that automatically controls the 
solution process also enables intelligent reanaiysls to incorporate sub- 
structure modifications. As input data describing the modifications are 
processed, substructure results dependent upon the modifications are in- 



17 


validated. For example, changing a load case definition invalidates ex- 
isting displacements for that load case but not the structure stiffness 
matrix. Prior to reanalysis, the system traverses the hierarchy to 
determine the type and order of computations required. During the 
traversal, dependent results for other substructures are destroyed and 
tagged for recomputation. When a substructure at level "I” Is marked 
for computation, dependent substructures at level "1-1" are also marked 
for computation. In this manner, the effects of substructure modifica- 
tions are automatically propagated upward through the hierarchy. POLO- 
FINITE equation solvers are based on the hypermatrix techniques first 
Introduced In the ASKA system. An efficient "partial decomposition” al- 
gorithm is utilized to condense substructure stiffness matrices. 

The most recently publicized substructure system, MISA [2.7], was 
developed by the Japanese ship building Industry. MISA Incorporates 
several unique concepts. Wavefront, rather than variable bandwidth, 
equation solvers are used to condense substructure stiffnesses although 
no computational advantage Is claimed. The software logic deduces the 
substructure hierarchy from parent-child relations Input by the user. A 
"copy" function facilitates repeated use of previously defined substruc- 
tures. MISA does not support mixed substructures and finite elements at 
the same level. It currently analyses linear structures for static and 
thermal loads, and performs steady-state heat conduction analyses. 

Dodds and Lopez C2.1H extended the POLO-FINITE system to support 
multilevel substructured models for static, nonlinear analysis. The 
analyst defines linear regions that are substructured and condensed to 
form effectively elastic supports surrounding the nonlinear, highest 



level structure. The reduced size of the nonlinear structure analyzed 
with the iterative technique- frequently yields significant cost savings. 
Currently, the nonlinear region must be defined prior to beginning the 
analysis. This is a major disadvantage of the approach. For the ini- 
tial loading levels, this proves inefficient as a much smaller nonlinear 
•region develops than is actually declared. To improve the situation. It 
is necessary to consider substructures that can be made nonlinear and 
brought into the iterative solution as the loading levels increase. 
This requires yet another level of sophistication in the control logic 
and has not been attempted. 

2.4 User- Software I nterf ace 

The user-software Interface for general purpose systems must 
provide sufficient flexibility to Invoke the options, and yet It should 
not discourage the Infrequent user with unnecessary details. Substruc- 
turing further complicates this Interface with the introduction of more 
elaborate topology, geometry, and computational algorithms. The 
developers of most software to support substructuring have not made the 
interface particularly simple for the user. These systems usually In- 
corporate substructuring as an extension of the original software 
design. Substructuring Is viewed as a last resort to obtain a solution, 
rather than as a natural modeling approach. Users are required to 
"program" the substructured solution by directly or Indirectly Invoking 
computational modules and by manipulating disk files that store sub- 
structure data. Consequently, only the most experienced users attempt 
substructured solutions; casual users are told to avoid substructuri ng. 
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The process of defining a substructured model and conducting the 
analysis consists of five logical steps. These are: 


1. The definition of each independent substructure (elements, 
constraints, and loadings). 

2. The elimination of substructure "slave" (condensed) nodal 
freedoms. "Master" nodal freedoms remain after condensation. 

3. The connection of individually defined substructures In a 
topologic and geometric sense. This requires matching of sub- 
structure boundaries to Insure displacement compatt Di I ity and 
may also include coordinate transformations of substructure 
matrices from their local system to a common global reference 
frame. Additional boundary constraints may also be required. 

4. The creation of a loading set hierarchy that parallels the 
substructure hierarchy. Equivalent nodal loads computed for 
the loading cases defined in (1) for the Individual substruc- 
tures are also reduced through condensation. Loads reduced to 
the master nodes serve to define a hierarchy of loading sets 
on the complete model. For example, the reduced loads may be 
applied to selected copies of a substructure at the next 
higher level to create a destred pattern of loading. The 
loads may also be carried up through the hierarchy with scalar 
multipliers and possibly combined with other reduced load 
cases to form new loading cases on a higher level substruc- 
ture. 

5. Requests for computation and output. The complexity of com- 
putational requests depends on the level of procedural detail 
required by the software to effect the solution. Output re- 
quests become complex when the user designates the hierarchy 
level for which results are desired, for which loading cases, 
elements, nodes, coordinate reference frames, etc. 


The common technique to approach these problems Involves the use of mul- 
tiple programs or a single program executed many times. The programs 
communicate through data stored on disk files (and/or tapes). The 
analyst is responsible for coordinating the program executions to 
produce the desired results. A typical analysis with one level of sub- 
structures might require that the following tasks be performed by the 
analyst: 
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1. Elements, nodal coordinates, topology, constraints and loads 
are defined to a finite element program for each individual 
substructure. The analyst provides special instructions in- 
dicating a condensation is desired and the list of master (or 
slave) DOF. Disk files are attached to the program onto which 
the reduced stiffness matrices and load vectors are written. 
Usually each substructure must be processed in a separate 
program execution and with unique data files. 

2. A substructure processing program is executed with all sub- 
structure disk files attached. This program is provided with 
disk file numbers for the substructures and the order In 
which they are to be processed. The connectivity relations 
between substructure master DOF and the global numbering 
system are also specified. Additional complexities arise when 
substructure reference axes are not all parallel. Since only 
substructures can be processed, it Is not possible with such a 
scheme to mix substructures and simple finite elements at the 
same level, which is Inconvenient In nonlinear analysis. The 
substructure processing program assembles the reduced stiff- 
nesses and loads to form a final set of equilibrium equations 
for solution, then computes the displacements for the highest 
level. Master DOF displacements for user selected substruc- 
tures are extracted and written onto another set of disk 
files. 

3. Finally, the finite element program is again invoked with 
proper disk files attached that contain displacements for the 
substructure master nodes. Backcondensation procedures may 
then be performed to recover slave displacements and element 
stra Ins-stresses. This process must be repeated for each copy 
of the substructure since the computed displacements, unlike 
the stiffness matrix, are unique. 


Each of the above steps Involve considerable manipulation of disk files 
and several program executions by the user. The manpower costs can 
easily approach those for model generation. The process is almost in- 
tractable when more than one level of substructuring Is used (if per- 
mitted at all by the software). Nonlinear problems require looping over 
the second and third operations for each load Increment. To place 
numbers on the amount of user effort Involved with this approach, a 
linear model containing four substructures at the first level was 
analyzed with one of the most successful commercial finite element 
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programs. Following the recommended procedures in the user documenta- 
tion, nine separate computer Jobs were required to obtain a solution. 
Absolutely no automatic solution techniques are available in the 
program. Changes to one of the substructures required a complete 
reanalysts of the entire structure due to the arrangement of data on 
tape storage. 

Demands placed on the user are greatly simplified when the software 
is designed to support substructured analysis In an integrated fashion. 
Standard analyses become the simplest default procedure in such a 
system. The following example, analyzed with the POLO-FINITE system. 
Illustrates the degree of simplification possible with a user-oriented 
approach to substructured analyses. This example Is not intended to 
demonstrate savings of computational effort; the computational advan- 
tages are demonstrated In more complex examples at the end of this 
chapter. 

The structure is a simple two span, planar truss as shown in 
Fig. 2.1. The generally non-symmetri c loads requires that a full model 
(both spans) be analyzed. Components of the substructured model are 
shown In Fig. 2.2, with names assigned to each component for identifica- 
tion In the model description. Figure 2.3 illustrates the substructure 
hierarchy In tree form for this simple example. One span is defined and 
condensed to the four nodes necessary for connection with the adjacent 
substructure. The bridge Is defined using two copies of the condensed 
substructure and three rod elements to complete the model. 
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The lowest level (substructure, SPAN, consists of 8 nodes (each 
with 2 DOF) and 13 rod elements. Input data describing structure SPAN 
to FINITE: are listed in Fig. 2.4. Element types, properties, topology, 
arid nodal coordinates are first specified. The problem oriented 
language (POL) input eases data entry by eliminating column and command 
.order restrictions. No natural boundary conditions occur at the nodes 
eliminated by condensation. Constraints are thus omitted on this lowest 
level substructure. Three independent loading cases are defined to act 
on SPAN. These represent a uniform load applied over the bottom chord, 
a uniform load on the top chord, and a lateral load acting on the top 
chord. The magnitude of each loading is unity to simplify the defini- 
tion of actual loading magnitudes in the higher level structure. 

Structure SPAN_C0N is defined as the statically condensed version 
of SPAN. Figure 2.5 lists the Input data describing this structure. 
Nodes 1, 3, 7, and 8 are retained after condensation. SPAN_C0N in this 
example corresponds to a "super-element" In the terminology used by some 
researchers. Analysts explicitly Introduce condensed substructures Into 
the hierarchy through intermediate structures such as SPAN_C0N. Struc- 
ture SPAM is referred to as the "child" of the "parent" structure, 
SPAN_CON, which resides at the next higher level in the structure tree. 
This technique has proven to be a natural means of Incorporating the 
condensed version of the substructure into the hierarchy. It eliminates 
confusion on the analyst’s part and maintains a consistent definition of 
structures in the database. Some structures are simply tagged as "con- 
densed" which serves to control execution of the processors. The in- 
cidences specified for element one of SPAN__C0N designate the nodes 
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C 

C PROBLEM UNITS ARE KIPS, FEET 

C 

STRUCTURE SPAN 

NUMBER OF NODES 8 ELEMENTS 1.3 

ELEMENTS ALL TYPE ROD E 4.32E06 AX 0.0347 

COORDINATES 

X Y 

10 0 

2 20 0 

3 20 20 

4 40 0 

5 40 20 

6 60 0 

7 60 20 

8 80 0 

INCIDENCES 

1 1 3 

2 2 3 

3 3 4 


LOADING UNITJTOP 
NODAL LOADS 

3 7 FORCE Y -10 
5 FORCE Y -20 
LOADING UNIT_B0TT0M 
NODAL LOADS 

246 FORCE Y -20 
LOADING UNIT_SWAY 
NODAL LOADS 
3 FORCE X 1.0 


Figure 2.4 — POLO-FINITE Input Data for Structure SPAN 
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STRUCTURE SPAN_CON 
NUMBER OF NODES 4 ELEMENTS 1 
ELEMENT 1 TYPE SPAN CONDENSED 
INCIDENCES 

1 1378 $ BECOME NODES 1-4 

LOADING UNIT_TOP 
EXTERNAL ELEMENT LOADS 
1 UNIT_TOP 

LOADING UNIT_BOTTOM 
EXTERNAL ELEMENT LOADS 
1 UNITJOTTOM 

LOADING UNITJSWAY 
EXTERNAL ELEMENT LOADS 
1 UNIT SWAY 


Figure 2.5 — POLO-FINITE Input Data for Structure SPAN_CON 
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retained from substructure SPAN. FINITE currently requires that all DOF 
at a node be either eliminated or retained during the condensation 
process. Loading cases are carried up through the structural hierarchy 
with a loading type designated EXTERNAL ELEMENT LOADS. An EXTERNAL 
loading specifies the names of loading cases on the child substructure. 
The loads are condensed and placed on the parent substructure nodes 
under the declared load case. There is a one-to-one correspondence 
between loading case names for the parent and child substructures in 
this example, but this Is not required. Only loading case names wtthin 
a structure must be unique. Any number of loading cases, with optional 
scalar multipliers, may be selected from the child substructure to con- 
struct loadings on the parent substructure. After the stiffness and 
loadings on SPAN are condensed during solution, SPAN_C0N has the major 
characteristics of any other structure; It has a stiffness matrix and 
loading cases stored In a standard format. 

Structure BRIDGE Is modeled from two copies of SPAN_C0N with three 
additional bar elements to complete the model. Figure 2.6 lists the in- 
put data describing structure BRIDGE. Copies of SPAN__C0N (elements 1 
and 2) are placed in BRIDGE with the same orientation relative to the 
coordinate system in which they are defined. In this model, substruc- 
ture reference axes X -Y and the highest level structure axes X a -Y are 
parallel. In other cases, a rotational transformation of substructure 
stiffness matrices and nodal loads may be required for proper geometric 
alignment. Coordinates are required for nodes of the three additional 
bar elements. Coordinate values are used to determine element size and 
orientation, thus the origin location is immaterial. For linear 
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STRUCTURE BRIDGE 

NUMBER OF NODES 8 ELEMENTS 5 
ELEMENTS 

1 2 TYPE SPANjCON ROTATION SUPPRESSED 
3-5 TYPE ROD E 4.32E06 AX 0.0694 

COORDINATES 

X Y 

2 0 0 

5 -20 20 

6 0 20 

7 20 20 

INCIDENCES 

1 14 5 2 

2 2 7 8 3 

3 5 6 

4 6 7 

5 2 6 

CONSTRAINTS 

1-3 V = 0.0 

1 U = 0.0 

LOADING FULL_TOP 
EXTERNAL ELEMENT LOADS 
1 2 UNITJOP 2 
NODAL LOADS 

5 7 FORCE Y -20 

6 FORCE Y -40 
LOADING LEFTJOTTOM 

EXTERNAL ELEMENT LOADS 
1 UNIT_BOTTOM 


COMPUTE DISPLACEMENTS FOR STRUCTURE BRIDGE 
[ output requests ] 

STOP 


Figure 2.6 — POLO-FINITE Input Data for Structure BRIDGE 
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analysis, mixing of substructures and elements Is more a convenience for 
model definition than a computational necessity. Constraints imposed on 
the nodes of structure BRIDGE model the simple support boundary condi- 
tions. EXTERNAL loads are again used to apply loadings from SPAN_C0N to 
BRIDGE. Load cases on BRIDGE are defined by selecting external loads on 
elements 1 and 2, Nodal loads, standard element loads (e.g., a 
distributed force on an element), and external element loads may be com- 
bined in any manner to define the loading cases on BRIDGE. 

A request for analysis has the simple form COMPUTE DISPLACEMENTS 
... as shown in Pig. 2.6. FINITE processors traverse the hierarchy to 
automatically determine the order of processing required and to check 
for topological consistency. The solution then proceeds to completion 
without user intervention. A solution In this context Implies the conr- 
putation of displacements for structure BRIDGE. The hierarchy traversal 
performed for each COMPUTE DISPLACEMENTS command insures that only 
needed quantities are actually generated. 

Output requests provide the capability to designate the unique 
copy of a substructure for which results are desired in addition to 
other Information such as loading cases, element and node numbers, and 
coordinate systems for stresses- stra I ns. The "substructure stack" is a 
part of the structure name In the OUTPUT command. Thus, to obtain 
stresses in all bar elements for the right span of BRIDGE for all 
loading cases, the following command Is sufficient; 
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OUTPUT STRESSES FOR STRUCTURE BRIDGE/2/1 

In which the 2/1 designates subelement 2 of structure BRIDGE, which is a 
structure named SPAN_C0N, then subelement 1 of SPAN_C0N, which Is a 
structure named SPAN. The list of elements and/or nodes (in this case 
an implied "all" elements) refers to the final structure listed in the 
stack. FINITE processors examine the substructure stack and 
automatically determine the type and order of backcondensation processes 
required to satisfy the request. 

The OUTPUT command provides considerable flexibility for requesting 
substructure results. The above command represents an extreme case In 
which the stack points to a unique occurence of a substructure In the 
hierarchy. Alternatively, the command 

OUTPUT STRESSES FOR STRUCTURE BRIDGE 

requests the prtnting of results for all elements in structure BRIDGE. 
In this case, FINITE processors automatically traverse the complete 
hierarchy below structure BRIDGE, recovering displacements, strains, and 
stresses at every level. The output processor traverses the same stack 
to print the results from the top (BRIDGE) down. 

The major advantages of the above approach to substructured model 
definition and solution are now apparent. The complete problem defini- 
tion and solution Is accomplished In a single computer run with no user 
Interaction concerning the placement of substructure Information onto 
data files. In effect, the analyst has access to a structural model 
database which can be defined and modified at a very high (logical) 
level. Substructure model description is a completely natural extension 
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of standard model description procedures. When implemented in this 
manner, computation requests are identical with and without substruc- 
turing. The traversal processes are performed automatically regardless 
of the model complexity. Output commands are quite flexible in that 
results for a wide or limited range of substructures may be requested. 
The POL input, while not required, considerably simplifies data entry 
and provides ready documentation for the model. It Is mandatory for in- 
teractive processing. 

FINITE currently requires that the structural hierarchy be defined 
in an inverted order. For example, structure SPAN must exist in the 
database at the time structure SPAN_C0N Is defined. Similarly, struc- 
ture BRIDGE cannot be defined unless structure SPAN_C0N exists. The 
structural tree must be defined from the bottom upj however, each branch 
need not be completely defined before beginning another branch. This 
restriction is sometimes inconvenient In that descriptions of a tree 
from the top down may be more natural. 

2 °5 Computational and Software Issues 

Multilevel substructuring creates a number of computational and 
software problems not encountered In conventional finite element 
systems. This section provides a brief survey of the major com- 
putational and software issues. The algorithmic details of stiffness 
condensation, load reduction, and displacement recovery for an Isolated 
substructure are first reviewed. The logical control techniques and 
data structures required to automatically process complex substructured 
models are then examined. These non-numericai aspects of the software 
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determine, to a large extent, the eventual generality and user- 
acceptance of the system. The final topic discussed In this section 
concerns data structures and algorithms for the solution of very large 
sets of linear equations. The less widely known hypermatrix procedure 
Is described. Hypermatrix computations offer a number of advantages 
compared to other techniques, including skyline and wavefront, espe- 
cially when the implications of a paged virtual memory system and 
parallel-pipeline hardware are considered. Hypermatrix techniques are 
also advocated later in this report for etgenproblem solution of very 
large systems. 

2.5.1 Substructure Reduction 

The three major computational activities associated with processing 
an Isolated substructure include: 

1. Condensation of the substructure stiffness matrices to 
eliminate the slave nodal freedoms? 

2. Condensation of the equivalent nodal loads applied to the 
slave nodes; 

3. Recovery of slave node displacements once master nodal 
displacements are known from solutions of higher level struc- 
tures. 

The classical equations for these operations were first presented by 
Przelmlenlckl C2.15] and are included here for completeness. The formal 
procedures are Inefficient and seldom followed. Two additional conden- 
sation procedures, referred to as coordinate transformation and partial 
decomposition, have better efficiency and are discussed In detail. 
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Equilibrium equations for an Isolated substructure are first par- 
titioned into two sets corresponding to the slave and master nodes as 


,ss 


,ms 


,sm 


,mm 


m 


( 2 . 1 ) 


where the superscripts m and s denote master and slave nodes respec- 
tively. The number of slave DOF is designated by "p", the number of 
master DOF by "q”. Constraints forcing a slave displacement to have a 
prescribed value and multi- point constraints between slave displacements 
are assumed to have been imposed by modification of the coefficient 
matrix and load vector(s). The solution of Eq. (2.1) in partitioned 
form follows the standard procedures summarized below: 


{u 5 } = [K SS r {p S > - CK SS l" 1 [K Sm ]{u m }, 


( 2 . 2 ) 


([K mm : - [K mS ][K SS ]" 1 [K Sm ]){u m } = {P Qff }, 


(2.3) 


{p eff } = {p m } - CK ms :cK ss ] ‘{p 5 } 


(2.4) 
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These equations reveal the form of the condensed stiffness matrix, 
Eq. (2.3), the reduced load vector, Eq. (2.4), and the procedure to 
recover slave node displacements from known master node values, 

S S I 

Eq. (2.2). Their inefficiency derives from the computation of [K J . 
The computational penalty Increases dramatically when [K SS H has a narrow 
•bandwidth but a fully populated Inverse (the most common case). Opera- 
tion counts for the computations are given later In this section for 
comparing the reduction methods. 

In the second method for condensation, the slave node displacements 
are related to the master node displacements through a coordinate trans- 
formation matrix, D"3» such that 

{u S } = CT]{u m }. (2.5) 

Each column of CT] contains the displacements of the slave nodes for a 
unit value of one master node displacement component, all other master 
DOF displacements being held zero. Because these displacements 
represent deformed substructure shapes that are analgous to mode shapes 
In dynamics, they are often referred to as "static constraint modes". 
The matrix Ct3 Is evaluated by substituting Eq. (2.5) Into the first row 
partition of Eq. (2.1) which yields. In the absence of external loading, 

[K SS 3CTD{u m } + CK Sm :{u m } = (0). (2.6) 

After eliminating (u m }frorn both sides of this equation, the product of 
[K s tl with each column of [TH equals the corresponding column of -CK sm Il. 
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Equation solving techniques are thus applicable to compute the columns 
of D"]. The matrix C« ss ] is triangulated only once using the front, 
Cholesk! or some other decomposition scheme. The procedure may be rela- 
tively efficient since [K ss ] is often banded (the natural ordering of 
slave DOF Is not altered in forming these equations). The columns of 
CT] are obtained by successive forward and backward substitution over 
the columns of -CK sri ]. An expression for the condensed stiffness matrix 
Is obtained by equating the strain energy of the substructure with and 
without the coordinate transformation. Thus, 

[K] = CT] T CK SS ][T] + [Tl T [K Sm ] + [K mS ][T] + CK™]. (2.7) 

Using the symmetry of off-diagonal submatrices, CK sm ] and C^ s ]» and the 
symbolic form of CT], the right side of Eq. (2.7) can be expanded and 
simplified. The form of the condensed stiffness for computation becomes 

CK] = [K mm ] + [K mS ][T]. (2.8) 

Eq. (2.8) Is often written In an unexpanded form that Includes an iden- 
tity matrix in the definition of CT]. The condensed load vector for the 
master DOF Is related to the slave DOF load vector by [T] In a similar 
manner using contragredience 


{p} = {P m } + Ct] t {p s } . 


(2.9) 



37 


A third technique for substructure condensation, which has been 
widely adopted, employs a "partial decomposition” of the stiffness 
matrix. Wilson [2.193 provides a detailed explanation of the com- 
putational procedures. Equations for an isolated substructure are as- 
sembled and partitioned into slave and master nodal DOF as in Eq. (2.1). 
•Gauss or the more efficient Choleskt decomposition is applied to com- 
pletely eliminate the first "p” rows correspond! ng to the slave DOF. 
Row-wise storage, and decomposition of the lower triangle Is assumed In 
this discussion. Similarly, the ma ster- slave coupling terms of CkT 1 S I] 
are reduced following standard procedures for off-diagonal terms. The 
Choleski reduction formulas applicable to the slave DOF are 

1-1 


. . = ( k . . - Z I ?. ) 1 /2 

m u k=1 ik 


i < P 


i J 


j-1 

k. . - Z 
' J 


k = i ' k J’ k • i, .< P 

j < i 


( 2 . 10 ) 


I . . 

J J 


Lower limits of k-1 on the summations imply a fully populated coef- 
ficient matrix. Extension to accomodate a variable bandwidth is 
straightforward. 

A partial decomposition Is then performed on the remaining [K mrr I] 
submatrix of master DOF coefficients to eliminate the coupling effect of 
slave DOF In the matrix [K mS 3. This Is accomplished by terminating the 
normal summations in Eq. (2.10) at column ”p" for the terms !,J > p. 
The master DOF stiffness terms are modified by this process to reflect 
the original coupling between slave-slave DOF and slave-master DOF en- 
tirely within the master-master DOF terms. The modified submatrix [K mrt !l 
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is the desired condensed stiffness matrix for the substructure. 

Final elimination of the master DOF occurs during solution ot the 
highest level structure. The basis for partial decomposition is that 
slave nodes within a substructure do not have a topological connection 
with nodes elsewhere in the hierarchy. Thus, standard decomposition of 
the complete structure stiffness, generated without substructuring, 
simply skips the summations that involve zero coupling terms between 
slave DOF in the equivalently substructured model. 

Condensation of substructure loading vectors Is accomplished 
through a forward reduction using the partially triangulated stiffness. 
Summations for master DOF terms are again terminated at column "p". 
Condensed loads for the master DOF reside In the last "q" rows of the 
load vectors. The first "p" rows contain partial displacements of slave 
DOF reaay for backsubstitutlon. These are termed the '‘partial slave 
displacements." 

Recovery of slave DOF displacements Is performed by completing the 
backward substitution over the first "p" rows once master DOF displace- 
ments are available for the next higher structure. Before the back sub- 
stitution is begun, the partial slave displacements generated and saved 
during load condensation are placed In the first "p" rows ot the load 
vector. Although this appears to be a trivial task, ft becomes a very 
complex logic problem when users are permitted to define load cases on 
higher level structures as combinations of condensed loading cases. In 
effect, users may implicitly define new load cases on the lower level 
substructure for which partial slave displacements are not computed 
during load condensation. The displacement recovery logic must 
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traverse the hierarchy again from the top, down to the level being 
processed to generate the loading combinations defined Implicitly. Only 
then can the correct partial slave displacements be computed by summing 
the separate load case values with the Implicitly defined multipliers. 

Computational efficiency is always of concern when large matrices 
.are manipulated. Experience with a large number of analyses has shown 
that generation of the condensed stiffness matrix requires nearly all of 
the computational effort associated with processing an individual sub- 
structure. Load case condensations and slave displacement recovery com- 
bined seldom require 10$ of the effort for sttffness condensation. Each 
of the algorithms described above Involve extensive summations to 
generate* resultant terms. The number of multiplications performed In 
evaluating these summations provides a measure of their relative ef- 
ficiency. Operation counts were developed for the stiffness condensa- 
tion phase of the three algorithms based on the following assumptions: 

s s 

1. Matrix [K J Is symmetric, has "p” rows, and an average half 
bandwidth of "r”; 

2. Matrix [K ms 3 = QK sn l! T and Is fully populated as a consequence 
of DOF reordering necessary to partition slave and master DOF 
for condensation; 

3. Matrix CK mri> ] Is fully populated. 

In the best possible case matrix CK ms ] Is upper triangular; the 

fully populated worst case is assumed here. Matrix QK 2 3 * S ^] will nearly 
always be banded. With these assumptions. It Is a simple task to 
estimate the number of multiplications required for each algorithm. The 
resuits are tabulated below: 
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1. Classical Inversion: p(0.5r 2 + pr + pq + q 2 ) 

? 7 

2. Coordi nate Transformation: p(0.5r + qr + pq .+ q ) 

? 7 

3. Partial Decomposition: 0.5p(r + pq + q ) 

To illustrate the effort required for each algorithm, consider a sub- 
structure for which p=500, q=50, and r=100; only 10$ of the nodal DOF 
are retained after condensation. The operation counts for each al- 
gorithm are listed below: 

1. Classical Inversion: 4.125 X 10 7 

2. Coordinate Transformation: 1.875 X 10 7 

3. Partial Decomposition: 0.940 X 10 7 

The results demonstrate clearly the Inefficiency of the classical 
inversion algorithm employed In early substructure software. The 
coordinate transformation algorithm, which is explicitly required for 
dynamic reduction, and the partial decomposition algorithm are con- 
siderably more efficient. The effort for coordinate transformation 
rapidly approaches that for classical inversion when "q" nears the value 
for "p", that Is, when a larger percentage of substructure DOF are 
retained after condensation. The results also demonstrate that the par- 
tial decomposition procedure should always be used for static analysis. 
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2 .5 .2 Logical Control and Data Structures 

The numerical algorithms described in the previous section are ap- 
plicable for an Isolated substructure. Before the computations for a 
substructure may begin, the software logic must determine the correct 
order In which to process the substructures. The proper ordering 
depends on the type of operations to be performed (such as stiffness as- 
sembly or displacement recovery) and the topologic relationships between 
substructures, as represented by the hierarchical tree. This Is not 
trivial task for models that consist of forty or fifty substructures 
distributed through five and six levels In a hierarchy. Finite element 
models with such characteristics are quite common for the analysis of 
large, geometrically complex structures. This section first outlines 
the familiar techniques used In standard analysis programs. An overview 
is then given of the techniques required to support general substruc- 
tured analysis as Illustrated in the POLO-FlNlTE example of Section 2.4. 

Programs that support only standard, static analysis have compara- 
tively simple control and data structure problems. A solution is accom- 
plished by Initiating a set of processing modules in a completely 
predetermined order that does not vary from one analysis to another. 
The typical sequence is: (1) read the input data, (2) compute the ele- 
ment stiffnesses, (3) assemble the structure stiffness and triangulate, 
(4) generate the loading vectors and perform a load pass to obtain 
displacements, (5) compute element strains and stresses, and (6) output 
results. Extensions to accomodate new loading cases are quite simple. 
Modification of the structural model necessitates a completely new solu- 
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tion. 

The data Is s+orea In a sequential manner that reflects the simple 
access mechanism required In the processing modules. A few sequential 
disk files suffice for most programs. Random files are occasionally 
adopted to facilitate stiffness assembly and trlangulatlon for large 
.structures processed with out-of-core techniques. 

In contrast, the automatic solution of a general substructured 
model requires a dynamic control capability. Each model solved with a 
general system may potentially require a unique order of module execu- 
tion. The order of execution cannot be "pre- programmed" In the software 
as It Is for standard analysis. Rather, the software logic must use a 
description of the substructure hierarchy to determine the flow of ex- 
ecution for each particular analysis (or reanalysis). Dynamic control 
logic easily accomodates modifications to substructures that alter the 
flow of execution during reanalysis. More sophisticated data structures 
and access schemes are needed to support the dynamic nature of the solu- 
tion procedure. The requirement for equal access to the data for any 
substructure eliminates consideration of sequential file storage (unless 
literally hundreds of files are available). The topologic dependencies 
between substructures suggest the use of a hierarchical data organiza- 
tion that parallels the natural substructure hierarchy. Formal data 
base management techniques are necessary to define, maintain, and access 
the data structures. With this approach it becomes feasible to ef- 
ficiently maintain data for any substructured model, regardless of size 
and complexity, within a single disk file. 
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For purposes of discussion, the substructured solution is con- 
sidered to have two computational phases. These are: 1) processing the 
model through the computation of displacements for the highest level 
structure, and 2) recovery of substructure displacements and element 
strai ns-stresses. The separation Into these two phases follows from the 
data structures and processing logic naturally suited for each task. 
Output generation is also Important but does not impact the com- 
putational processes or data organization. The first phase involves 
stiffness and load vector assembly, and requires that substructures be 
processed upward from the bottom of the hierarchy. Data structures must 
support the repeated use of a .substructure stiffness matrix and loads to 
form similar matrices for higher level structures. In the second com- 
putational phase, processing of the hierarchy occurs from the top-down 
and then only along user designated paths. Displacements, strains, and 
stresses are not normal ly recovered for a 1 1 substructures. Data struc- 
tures must reflect the uniqueness of displacements, strains, and 
stresses for each occurence of a substructure in the hierarchy. 

Stiffness matrix and load vector assembly proceed upward from the 
lowest levels of the hierarchy. Assemoly of a structure at level "I" 
cannot begin until all dependent substructures at level ”1-1" have been 
completed. Control of the assembly process through a stack driven 
traversal of the hierarchy appears, at first, to be a natural approach. 
However, a topologlc sort of the hierarchy (performed using a stack) to 
determine the processing order before computations begin is far simpler. 
The topologlc sort provides the stiffness module with a properly ordered 
list of substructures to process. Figure 2.7 presents a block flow 
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Figure 2.7 — Topoiogic Sort to Order Substructures for Assembly 
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diagram for the topologic sort. A substructure appears only once in the 
sorted list regardless of its number of occurrences in the hierarchy. 
During reanalysis, the topologic sort simply omits ail unmodified sub- 
structures from the list. Logical control within the assembly module 
becomes quite simple given the sorted list of substructures to process. 
Figure 2.8 illustrates the overall assembly logic, with the topologic 
sort shown as the first operation. The first structure in the list is 
extracted and its stiffness assembled, followed by the second structure 
in the list, etc. Only structures appear In the sorted list. The as- 
sembly module generates matrices for finite elements which occur in the 
structure being assembled. 

When a ’’condensed” substructure is extracted from the sorted list, 
the assembler will have already computed the "uncondensed” stiffness 
matrix since the single child of the condensed substructure occurs one 
level lower In the hierarchy. The assembler suspends execution and in- 
vokes the equation solving module, passing it the identifier for the 
substructure to be condensed as Indicated In Fig. 2.8. After the stiff- 
ness condensation is completed, the assembler resumes execution, 
retrieves the master node matrices, and saves them In data structures 
for the condensed structure. No traditional assembly operations are 
performed. Thereafter, the condensed substructure has a stiffness 
matrix with the same format as any other structure (or finite element). 
The assembler continues execution by selecting the next structure in the 
sorted list. Because the processing modules are dynamically invoked, 
condensation at any number of levels creates no logical control dif- 
ficulties with this procedure. 
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Figure 2.8 -- Multilevel Substructure Assembly Procedure 
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Figure 2.9 illustrates the essential features of a data structure 
for storage of the hierarchy description and stiffness matrices. The 
ELEMENTS table contains one column for every finite element and struc- 
ture declared by the user. Rows of this table describe attributes for 
the element or structure, for example, the number of nodes "NUMJ'IODES". 
A STRUCTURE table exists only for columns that contain a structure. The 
STRUCTURE table stores pointers to lower level tables that contain data 
applicable only to a structure, for example, subelement incidences, 
subelement orientations, nodal coordinates, and components. The COMPO- 
NENTS table has one entry for each subelement in the structure. The 
data value for each subelement is the column number, denoted ECOL on the 
figure. In the ELEMENTS table that defines the subelement. The subele- 
ment may be a structure or a finite element. Given the ECOL for a 
structure, the hierarchy from that level downward may be easily 
traversed. 

The stiffness matrix for a finite element and for a structure have 
Identical storage formats as shown in Fig. 2.9. A stiffness matrix con- 
sists of nodal "submatrices" which are stored In a row- wise format. The 
column of the first non-zero submatrix in a nodal row Is Indicated by 
the value of FIRST_C0L. For finite elements, this value Is always one. 
For structures, the value of FIRST_C0L reflects the half bandwidth for 
the nodal row. The lower triangle of the element stiffness Is stored. 
All submatrices from column FIRST_COL through the diagonal are generated 
even though there may be Intermediate null submatrices. 
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Assembly of the equivalent nodal loads for a multilevel substruc- 
ture model Is logically more complex than stiffness assembly. For ex- 
ample, when displacements for a load case defined on the highest level 
structure are requested, the equivalent loads processor must traverse 
the complete hierarchy while determining which substructures have load 
cases that contribute to the one specified on the highest level struc- 
ture. The procedure is even more interesting when analysts define load 
cases in terms of other load cases from one substructure level to the 
next. A more flexible data structure than a simple vector of structure 
identifiers Is necessary to accomodate the load case Identifiers. A 
multiple vector. Inverted list proves adequate. Condensation of the 
load vectors Is performed in the same manner as stiffness condensation — 
as an Interrupt In the normal procedure of transferring loads from one 
substructure level to another. 

Phase two computations recover substructure slave displacements and 
element stra I ns- stresses. The analyst generally specifies a path 
througn the hierarchy that Identifies the occurence of a substructure 
for which results are desired. The particular substructure may be at 
the lowest level of the hierarchy. In which case relatively few results 
are generated. Alternatively, it may be an Intermediate level substruc- 
ture, which can be used to Imply displacement recovery for that level 
structure and all lower level substructures of the same branch. 

Logical control using a stack technique is most convenient to drive 
substructure displacement recovery. Given the analyst supplied path and 
the COMPONENTS table shown In Fig. 2.9, the traversal procedure to reach 
the desired substructure (s) is quite simple. The back-condensation com- 
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pu+ations to retrieve slave displacements from the master values are 
logically treated as an Interrupt In processing down the hierarchy. In 
a manner paralleling the assembler condensation procedures, the 
displacement recovery module suspends execution, then Initiates the 
equation solving module to compute slave DOF displacements for the cur- 
rent substructure. On resuming execution, the displacement recovery 
module extracts the slave DOF displacements In the solver data struc- 
tures and reformats them to conform with standard displacement data 
structures. The next lower level In the substructure list may then be 
processed. Element strain-stress computations may be performed Im- 
mediately after slave DOF displacements are recovered or a separate 
module may be Invoked following the completion of all displacement 
recovery. A stack driven procedure also facilitates the computation or 
substructure stralns-stresses. 

The most Interesting aspect of phase two Involves the development 
of data structures for storage of substructure displacements, strains 
and stresses. Only one stiffness matrix exists In the database for a 
substructure regardless of Its number of occurrences, but the displace- 
ments, strains, and stresses are unique for each occurrence. The data 
structure shown In Fig. 2.9 for phase one processing requires extension 
to support substructure displacement recovery. One solution for this 
problem Is Illustrated In Fig. 2.10. The uniqueness of substructure 
results Is recognized only at the highest level structure. The LOADS 
table points to lower level tables that contain the definition and com- 
puted results for each load case. Displacements, strains, and stresses 
for the complete hierarchy for a load case are stored under the cor- 
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responding column of the LOADS table (defined for the highest level 
structure). Results are packed into sets of variable length vectors 
that reside under the DISPLACEMENTS, STRAINS, and STRESSES pointer 
tables. The allocation and retrieval of result vectors is accomplished 
by a two- level pointer scheme. For example, vector SPOINT shown in the 
•figure contains one entry for each subelement in the structure. For 
subelements that are finite elements, the "f” th entry in SPOINT 
defines the vector and position within the set of result vectors at 
which values for the element begin. For subelements that are substruc- 
tures, the SPOINT entry defines a relative shift in the vector numbers 
from the current level structure’s vector to the first relative vector 
for the specific occurrence of a substructure. SPOINT entries for the 
substructure refer to vector numbers that are relative to the first 
result vector for the substructure. Results for multiple occurrences of 
the same substructure thus appear in different stress-strain vectors. 
The absolute vector number for processing at any time is obtained by ac- 
cumulating relative shifts for each substructure present in the traver- 
sal stack. 

The most complicated aspect of this scheme Involves construction of 
the SPOINT vectors for each substructure. Note that the SPOINT vectors 
are Independent of any load cases. The mapping of substructure results 
onto vectors is identical for each load case. An identical scheme is 
used to store substructure nodal displacements, reactions, initial 
strains-stresses, and residual loads in nonlinear analysis. 
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2 „5 .3 Linear Equation Solving 

Efficient procedures for linear equation solving are essential In 
finite element analysis. Reduction of the symmetric stiffness matrix to 
triangular form often requires 50 percent of the computational effort 
for linear analysis and a slightly smaller percentage for nonlinear 
analyses., Many efficient algorithms and the details of their computer 
Implementation have been published. Meyer jl2 . 1 2 H has presented an ex- 
tensive review of the subject. 

The growing availability of computer hardware designed specifically 
for "number crunching" applications has spurred renewed Interest In 
equation solving procedures. The details of data storage and access 
mechanisms comprise the key factors In the ability to utilize advanced 
hardware. Advanced machine arch I tectures include parallel processors 
(the basis for most supercomputers) and virtual memory superminis. Some 
of the most recent supermini computers also have a pipeline design that 
greatly speeds up numerical operations compared to more conventional 
scalar processing. The pipeline concept uses a single processor but 
with separate prefetching of data to minimize processor wait times. It 
thus represents a intermediate design between a scalar processor and the 
parallel processing machine. The scalar design superminis can also be 
enhanced with the addition of an attached array processor. 

The basic requirements for any equation solving procedure can be 


summarized as: 
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1. The number of equations should be limited only by the amount 
of disk storage available; 

2. Use a reasonable amount of memory so as not to severely impact 
scheduling of the computer runs; 

3. Minimize the storage of zero coefficients and operations on 
them; 

4. Minimize data transfers to and from disk (I/O); 

5. Exploit capabilities of modern computer architecture. 

Substructure condensation by partial decomposition is essentially an 
equation solving procedure and therefore the same requirements listed 
above apply. However, stiffness condensation may adversely effect one 
decomposition procedure more than another. With band-based solvers, for 
example, reordering of the equations to facilitate condensation may 
produce a border-banded matrix of coefficients. The computation time 
required for partial decomposition is greatly increased compared to the 
time required for complete decomposition of the equations without 
reordering. 

This section briefly surveys three standard equation solving 
procedures, namely (1) Choleski with variable band (skyline) storage, 

(2) the front method of solution, and (3) Choleski with hypermatrix 
storage. Each of these Is a direct solution algorithm that involves a 
three step process: (1) reduction of the coefficient matrix to 
triangular form, (2) forward reduction of the load vector(s), and 

(3) recovery of displacements by a backward pass. The primary interests 
in the current discussion Include the Impact of reordering the equations 
prior to condensation and the prospects for adopting each procedure to 
new computer hardware. 
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Variable Bandw ! dth Procedure 

This tecnnique uses a CholeskI or Crout tr i angu I at ion algorithm 
with a compact storage scheme to accomodate wide variations in band- 
width. Jennings and Tuff [2.83 originally developed the procedure. 
Mondkar and Powell [2.13] later published the implementation details and 
results of timing studies. This procedure appears to have been widely 
adopted for "In-house" finite element systems because the data handling 
details are straightforward. 

Contiguous columns of the upper triangle are grouped together In 
blocks (Fig. 2.11). Ail coefficients In a block are transferred to and 
from disk In a single, logical operation. The available memory deter- 
mines the number of columns assigned to each block. The first non-zero 
row in each column is recorded during assembly. This Information Is 
simply appended after the last column of each block. All terms from the 
first non- zero row to the diagonal in a column are assumed to flll-tn 
during trtanguiation. 

Reduction of each term in column "I" requires that column "I 11 and 
at least one other column, "j", where J < I, reside In memory. Thus, at 
least two blocks of columns must fit simultaneously into the available 
memory. During forward load reduction, one coefficient block and one 
load block are required In memory (multiple load cases are handled very 
efficiently). The same memory requirements exist for the backward pass 
to recover displacements. 
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Over 80 percent of the effort in triangulating the matrix involves 
computing Inner products of two non-conti guous columns. This makes the 
algorithm very attractive for all types of advanced hardware. The 
number of terms In the inner products is usually large relative to the 
break even point on most hardware (some overhead is Involved In parallel 
processing to initially align data which requires that inner products 
have a minimum length to recover the overhead time — generally 5 to 10 
terms) . 

The variable band procedure has a number of drawbacks. It can be 
very I/O Inefficient during structure stiffness assembly if only one or 
two blocks fit Into memory. Elements are usually processed In sequen- 
tial order to assemble the structure stiffness. As a consequence, it Is 
possible for a structure to have a very narrow bandwidth and yet require 
considerable swapping of blocks during assembly. The classic example Is 
a narrow rectangular grid with nodes numbered in the short direction to 
minimize the bandwidth and elements numbered in the long direction 
(which does not affect the bandwidth). To assemble each sequence ot 
elements in the long direction, each coefficient block must be brought 
Into memory. Another drawback occurs when column heights vary sig- 
nificantly. Large numbers of inactive terms are transferred Into memory 
during triangulation as a consequence of the blocking (Fig. 2.11a). 
Lastly, the storage scheme Is biased towards accessing coefficients 
column-wise, which Is natural for triangulation and forward load reduc- 
tion. As noted above, each block Is transferred to memory only once 
during the forward pass. During the backward pass, however, row-wise 
access is the most natural. A block of columns may be transferred Into 
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memory many times depending on bandwidth variations. 

When executed on virtual memory hardware, the logical memory space 
(that declared in FORTRAN DIMENSION statements) for blocks may be very 
large, for example, sufficient to store seven or eight blocks. The 
program issues fewer logical I/O requests to transfer blocks than a con- 
ventional approach in which space for only two blocks is allocated. The 
virtual memory operating system performs additional I/O required forcing 
the program to execute In a given amount of real memory (working set 
size). This process Is normally handled through fragmentation of data 
arrays into pages which are transferred to and from memory as the array 
elements are referenced in the program. Considerable experimentation is 
necessary to balance the virtual paging and program block I/O, espe- 
cially on virtual machines that permit users to explicitly declare the 
working set size. 

Condensation causes difficulty only when the bandwidth increase due 
to reordering becomes excessive. Consider a structure which is numbered 
to minimize the average bandwidth for conventional analysis. In 
general, the reordering to accomodate condensation produces the border- 
banded matrix shown in Fig. 2.11b. This procedure requires that at 
least one complete column fit Into a block. For large 3-D structures, 
this requirement may lead to large increases in block sizes to ac- 
comodate the very long columns corresponding to master degrees of 
freedom. During the partial decomposition process, no computational 
penalty occurs until the equations with the large band in the last 
blocks are reached. Most of the triangulation time may be involved in 
eliminating these terms. However, advanced hardware can be exploited to 
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the fullest In processing the very long columns. 

The highly touted frontal procedure was first Introduced by Irons 
C2.6] and has since been extended to Include condensation. The frontal 
solver Is essentially Gauss elimination with extensive bookkeeping to 
minimize operations on zeroes. The structure stiffness matrix Is never 
explicitly assembled. Rather, element stiffnesses are brought sequen- 
tially Into memory, their terms added to the system matrix and then 
triangulated, all In one logical (but very complex) process. The 
ordering of equations In the system matrix Is determined by the element, 
rather than the node, numbering scheme. As elements are processed, 
newly appearing DOF are simply appended as new equations. The memory 
occupied by active equation coefflcents Is referred to as the "front.' 9 
The front storage space varies dynamically during solution as new equa- 
tions are added and old ones are completely eliminated (eliminated coef- 
ficients are transferred to disk). The front Is very similar to the 
"active triangle" concept In band-based solvers. 

The frontal procedure eliminates much of the CPU and I/O costs as- 
sociated with structure stiffness assembly. The trlangulatlon aspect Is 
no more efficient than the CholeskI process. There are two major disad- 
vantages of the frontal procedure for non-substrucfured analyses. 
First, the active front size becomes quite large for 3-D solid analyses. 
The bookkeeping logic Is very complex when the complete front fits Into 
memory; It appears Intractable when coupled with a "spill" algorithm to 
accomodate a large front that Is partially memory resident. Storage 
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space of 100-200K words for an active front Is common. This may cause 
scheduling problems on some machines and Impose absolute structure size 
limitations on others. The second difficulty with the frontal method is 
Its adaptation to advanced hardware. Uni Ike the Choleski procedure, 
inner products Involving hundreds of terms are not common In the frontal 
.solver. Operations are performed almost randomly over the active front. 
Parallel and pipeline hardware offer little advantage. 

The frontal method appears Ideally suited for virtual memory com- 
puters with scalar processors. The problem of handling a very large ac- 
tive front is relegated to the virtual memory operating system of the 
computer. An extremely large array space Is dimensioned for the front. 
The operating system pages the array segments as needed to maintain a 
predefined working set. The random accessing Into the front causes no 
penalty on a scalar (sequential) processor. 

The absence of an explicitly assembled structure stiffness Is a 
drawback of the frontal solver for substructured analyses. Substructure 
stiffness matrices (original and reduced forms) must always be available 
for use In defining other structural hierarchies. Stiffness condensa- 
tion also affects the efficiency of the frontal solver as It does the 
variable band solver. Consider the element grid shown In Fig. 2.12a. 
All nodes on the perimeter are retained after condensation. Conven- 
tional node numbering for a band solver produces a border-banded matrix. 
Condensation In the frontal solver Is accomplished by retaining all 
master nodes in the front until slave nodes are completely eliminated. 
For the mesh shown, the active front size Increases each time a top and 
bottom row element is processed. The front size reaches a maximum when 
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(a) Worst Numbering for Condensation with Front Procedure 
— Best Numbering for No Condensation 
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Figure 2.12 — Effects of Element Numbering on Condensation by Frontal Technique 
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element 88 is processed. Alternatively, elements can be numbered as 
shown In Fig. 2.12b for front solution. The front size remains constant 
while all slave DOF are completely eliminated. It then grows rapidly to 
the same maximum size for the previous numbering scheme as elements 
coupling the master and slave DOF are processed. The only advantage of 
the second case is that a maximum front size exists for less execution 
time. This could drastically reduce the paging rate on virtual memory 
machines. The same effect occurs for the variable band solver In that 
the active triangle reaches a maximum size when the master DOF are en- 
countered. 

Hypermatr ix Procedure 

The basic concept of hypermatrix storage is illustrated in 
Fig. 2.13a. A block of contiguous columns defined in the variable band 
procedure Is further partitioned row- wise to form rectangular "hyper- 
matrices". This storage format overcomes the major problems with 
variable band storage — excessive column heights and the transfer of 
unused terms during triangulation. 

Diagonal hypermatrices are always square; off-diagonal matrices are 
frequently rectangular. Hypermatr! ces with all zero terms are never 
created. It Is thus a simple matter to omit large numbers of zero 
operations with such a storage format. The sizes of hypermatrices can 
easily be adjusted to fit a particular amount of memory available for 
processing. The use of larger blocks (I.e., greater than 50 x 50) in- 
creases the number of zero terms picked up at the periphery of the band. 
The use of smaller blocks Improves the recognition of zero terms but In- 
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creases the data handling overhead. 

A formal data structure for hypermatrix storage is shown In 
Fig. 2.13b. Individual hypermatrices are stored as separate entities 
(for example, a logical disk file record). Pointer vectors that locate 
all hypermatrices within a column emlnate from the header table shown. 
■Pointers to zero hypermatrices are not maintained in the vectors. In- 
stead, the header table contains an offset locating the first non-zero 
hypermatrix in each column. 

Hypermatrix partitioning has been studied extensively by computer 
scientists [2.103 In connection with array operations on virtual memory 
(paging) hardware. It has been demonstrated that hypermatrix par- 
titioning requires the least working storage and Incurs the least page 
faults for the operations of multiplication and triangulation. The 
usual scalar formulations for matrix mul tipi Ication and triangulation 
(Gauss, Cnoleski, Crout) have very similar counterparts when cast In 
hypermatrix form. Choleski decomposition In hypermatrix form requires 
the same number of floating point operations as in scalar form. Both 
multiplication and triangulation require that only three hypermatrices 
reside simultaneously In memory. Therefore, large problems may be 
solved In a very small memory space, for example, 7500 words if 50 x 50 
blocks are used. Moreover, if memory for only three blocks is available 
it becomes a trivial problem to predict exactly the number of block 
transfers to and from disk. When space is available In memory for more 
than three blocks, it has been suggested [2.10] that a "least recently 
used” replacement algorithm efficiently utilizes the additional memory. 
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Fuchs and Schrem [2.43 devised a CholeskI triangula+ion scheme in 
hypermatrix form and implemented the procedure in the ASKA program 
[2.16j. The data structure in ASKA employs a two level matrix of 
pointers to the hypermatrices instead of the scheme shown in Fig. 2.13b. 
The POLO-FINITE system uses the hypermatrix procedure with the storage 
format shown in Fig. 2.13b. Both ASKA and POLO-FINITE use demand 
paging, virtual memory management systems for which hypermatrices are 
ideally suited. Hypermatrix sizes are usually defined such that one 
hypermatrix fits onto a page. A page corresponds to a single random 
disk file record. The advantage of allocating one hypermatrix per page 
Is that totally unbiased access to any block is obtained. Triangula- 
ti'on, forward pass, and backward pass operations thus have equal ef- 
ficiency with respect to data access. As previously noted, the variable 
band scheme suffers a heavy access penalty during the backward pass. 

As with the variable band procedure, the major effort in hyper- 
matrix triangulation involves Inner products of non-conti guous columns. 
Figure 2.14 illustrates this. However, in hypermatrix form, the inner 
product Is actually a sequence of matrix multiplications as indicated on 
the figure. In effect. Inner products over multiple scalar columns 
proceed simultaneously. Timing studies for large sets of equations have 
shown that the matrix multiplications usually represent 90 % of the 
triangulation effort. The data handling overhead seldom exceeds 5 % of 
the triangulation time. For example, after two 50 x 50 blocks are 
brought into memory, approximately 250,000 floating point mul tip! Ica- 
tions and additions (plus subscripting) are performed to complete the 
matrix mul tipi ication. 
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Hypermatrix tri angulation appears well suited for adaptation to 
parallel and pipeline hardware. The order of matrices in the multiply 
operation shown in Fig. 2.14 is particularly important. Because the 
first term is transposed, the matrix multiply Is simply inner products 
of columns In the two matrices, rather than the traditional row mul- 
tiplied into a column. Storage of the lower triangle in hypermatrix 
form eliminates this advantage of column-wise inner products. The 
equivalent matrix product requires inner products of two rows, rather 
than two columns. The choice of lower or upper triangle storage is Im- 
material for computations performed on scalar hardware. However, upper 

triangular storage is preferred for the most general case. 

The large bandwidth Increase that occurs with equation reordering 

prior to condensation presents no difficulties In the hypermatrix 
scheme. Large bandwidth fluctuations simply increase or decrease the 
number of non-zero hypermatrices in a column. The problem of fitting an 
entire column Into memory that occurs with the variable band method does 
not occur with hypermatrix storage. 

2.6 Examples of Substructured Analyses 

Two example analyses are described In this section to Illustrate 
typical computational savings achieved with substructured models. In 
the first example, a linear analysis Is performed for a portion of a Jet 
engine exhaust duct. Condensation is applied at two levels In the sub- 
structure hierarchy — at the lowest level to reduce the bandwidth and 
at the next higher level to utilize repeated components. The exhaust 
duct represents an excel lent test case for substructured dynamic 
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analysis with modal synthesis (discussed in Chapter 3) considering only 
linear response. Dynamic analyses using a full, unsubstructured model 
and some limited experimental data are available in the literature for 
comparison. 

The second example problem consists of a slightly curved thick 
shell constructed of impact resistant acrylic. A uniform pressure is 
applied over a very small area at the apex to simulate the Impact of a 
nondeformable object. The magnitude of the pressure Is Increased In 
each load increment but the loading area remains constant. Nonlinear 
response due to yielding of the material in the impact region is con- 
sidered. The problem typifies a large class of structures for which 
substructured models reduce the computation time. The region of non- 
linear response is easily estimated prior to the analysis. Standard and 
substructured models for this example are analyzed to provide data for 
comparisons of computational effort. This problem also provides an ex- 
cellent test case for substructured, nonlinear dynamic analysis. A 
transient analysis is required to predict dynamic response following the 
impact of a high velocity projectile. Even under such loading, the non- 
linear zone remains small relative to the overall structural dimensions. 
Condensation of the linear region should greatly reduce the com- 
putational effort for a transient analysis. For comparison, the com- 
putational effort within each time step of a transient analysis cor- 
responds to that for a load increment In static analysis. The results 
presented here for the static solution provide a basis to estimate com- 
putational savings for transient analysis. 



69 


Each structure has been analyzed with the POLO-FINITE system. Com- 
parisons between substructured and standard model solutions are made on 
the basis of CPU time and I/O. The POLO supervisor performs database 
and memory management functions for FINITE subsystems. POLO has exten- 
sive instrumentation that provides detailed summaries of CPU usage among 
the FINITE subsystems, for example: input, assembly, and triangulation. 
Within each subsystem, the CPU time expended on data management ac- 
tivities and on actual finite element computations is also available. 
CPU times are presented in non-dimensional form to eliminate dependen- 
cies on the processor execution speed. I/O activity is measured by the 
number of "page faults" executed by the POLO memory management system. 

A page fault represents the transfer of one 2500 word record from memory 
to disk followed by the transfer of another 2500 word vector to memory. 
For a given structural model, the number of page faults performed by 
POLO is independent of the computer hardware; It depends only on the 
dimensioned length of a data vector within POLO. Page faults are thus a 
very simple measure of I/O activity. Printer output and I/O transfers 
to sequential card image files are ignored. 

2.6.1 Lin ea r, Exarap.La 

Figure 2.15 shows a portion of a jet engine exhaust duct modeled 
for static stress analysis. The duct consists of a thin circular shroud 
connected to the central core with radial fins. Loadings of Interest 
include torsion, uniform external pressure, and nonuniform temperature 
distributions. For this analysis, the Inner edge of each fin is com- 
pletely fixed. 
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DIMENSIONS 
Shroud Fins 

R 3 170 mm L ■ 64 mm 
t = 1.5 mm t = 3.0 mm 

w = 81 mm w - 27 mm 



Figure 2.15 — Jet Engine Exhaust Duct for Example Linear Analysis 
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Difficulties in constructing a proper finite element model arise at 
each junction of a fin with the shroud. Accurate determination of the 
bending stresses in the junction vicinity requires a relatively fine 
mesh and the maintenance of element displacement compatibility. Curved 
shell elements adequately model major portions of the fins and shroud. 
The 8 node isoparametric shell element is used here. Each junction is 
modeled as a 3-D solid with shel l-to-sol Id, and sol I d-to-sol i d transi- 
tion elements employed to maintain displacement compatibility. This 
model provides realistic predictions of stress distributions at each 
junction without an undue increase In the number of nodal DOF. Figure 
2.16 Illustrates the four types of elements employed in the analysis. 

A 360 degree model Is required for general analysis. For the grid 
shown In Fig. 2.15, the full model has 408 elements and approximately 
9500 DOF. The generation of each element stiffness matrix requires 
numerical Integration. And while the large number of DOF does not 
present major problems for general purpose systems, the grid topology 
does introduce some severe computational penalties. The traditional 
node numbering scheme follows the narrow (axial) direction then the cir- 
cumferential direction. Each fin causes a large re-entrant area In the 
equation coefficients, but the most severe penalty arises because the 
first and last nodes are coupled. As a result the stiffness matrix is 
border- banded; the last 78 rows have a bandwidth of nearly 9500. This 
drastically increases the triangulation effort. Substructur i ng with 
condensation proves very computationally efficient for this structure as 
demonstrated below. 



ode Shell (QSHELL) 


b) Shell-to-Solid Transition 
(SHELL-TRANS) 



Solid (TS16! SOP) 

Figure 2.16 — Element 



d) 18 Node Solid Transition (TS18ISOP) 


Types for Exhaust Duct Model 
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Division of the full model into identical 90 degree substructures 
Is obvious. The selected partitions place the fin at the center of the 
substructure as shown in Fig. 2.17. The location of each type of finite 
element is also Indicated on this figure. Condensation of this sub- 
structure reduces the number of nodes to 26; 13 nodes across each end 
for connection with adjacent copies of the same substructure. Within 
the 90 degree substructure, the fin is also modeled as a substructure 
and condensed to the 13 nodes that connect to the shroud. Condensation 
of the fin eliminates the re-entrant area in the coefficients described 
above for the full model. Figure 2.18 shows the same 90 degree section 
model without substructuring. 

Considering the analysis of an Isolated 90 degree section, the fin 
condensation reduces CPU time by 14$ and I/O transfers by 21$ compared 
to an analysis without fin condensation. (Fin stiffness generation and 
condensation time are Included in the comparison.) However, condensation 
to eliminate all but the 26 boundary nodes increases the CPU by 53$ and 
the I/O transfers by 65$. These Increases reflect directly the penalty 
incurred by reordering the equations prior to condensation. Reordering 
for this substructure has the same effect on the bandwidth as closing 
the ring on the full model. The first and last nodes are effectively 
coupled by the reordering, but final 156 rows have a bandwidth of only 
2400 compared to 9500 for the final 78 rows of the 360 degree model. 
The large bandwidth difference produces the major savings for the full 
model generated with substructures. (Recall that the triangulation ef- 
fort Is proportional to bandwidth squared.) 
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13 Nodes on this Edge Retained 



S u bstructure Quarter 


Figure 2.17 — Substructured 90° Model for Exhaust Duct 
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Figure 2.18 — 90° Section of Exhaust Duct Without Substructur i ng 
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Using four properly oriented copies of the reduced 90 degree sub- 
structure, the 360 degree model was generated and solved. The final 
structure has only 312 nodal DOF with effectively a full bandwidth. 
(The condensed substructure stiffness Is fully populated.) Figure 2.19 
Illustrates the substructure hierarchy for the 360 degree model. Five 
levels of substructures are present In the model. 

Table 2.1 compares the CPU and I/O for the various solutions. The 
CPU and I/O required for complete solution of the substructured 360 
degree model are assigned values of 1.0. Relative times are given for 
other models. Over 90$ of the full model solution effort Is expended in 
the substructure reduction process. This is not too surprising given 
the very large DOF reduction and the very small number of DOF at the 
hignest level structure. It was not possible to analyze the unsubstruc- 
tured 360 degree model with POLO-FINITE given the current limitation on 
the number of nodes In a structure (the system limits the number of 
nodes in a single substructure to 833). However, accurate estimates of 
the CPU and I/O time are possible using the exact timings obtained for 
solution of the isolated 90 degree structure without fin condensation. 
The estimated CPU and I/O consider the factor of 4 In element stiffness 
generation effort, the factor of 4 In number of DOF, and and the 
quadratic Increase with bandwidth for the final 78 equations. The rela- 
tive savings with the substructured model are in the 13-14 range for CPU 
and 18-20 for I/O transfers. Most of the saving results from the band- 
width reduction noted above, but a minimum factor of 4 Is realized from 
the reduction In nodal DOF and element stiffness generation. 




Figure 2.19 — Structural Model Hierarchy for Exhaust Duct 



( 1 ) 

90° Standard 

0.69 

0.63 

(2) 

90° Substructured 

0.60 

0.54 

(3) 

90° Condensation 

0.92 

0.89 

(4) 

360° Substructured 

1 .0 

1 .00 

(5) 

360° Standard 

* 

13.6 

18.92 


* 

Estimated using CPU and I/O measures for the 90° 
standard model 


Table 2.1 — CPU Times for the Exhaust Duct Analysis 
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Figure 2.20 provides a detailed breakdown of CPU among the various 
POLO-FINITE subsystems for solution of the 360 degree substructured 
model. Mot surprisingly, the major effort is expended in the assembler 
and equation solver modules. It is Interesting to note that a rela- 
tively small percentage of the total job time (1455) is expended in data 
management activities. The percentage of data management time in the 
assembler is larger due to the very large number of small matrices that 
must be manipulated. In contrast, the equation solver accesses data in 

larger blocks (50x50) and performs a significant amount of numerical 
computation on them. Consequently, percentage of data management time 
Iri the solver is very small. 

2,6.2 Non I ? near ExampJja 

Figure 2.21 shows one quarter of a thick, shallow shell structure 
constructed of impact resistant acrylic. A monoton I ca I I y Increasing 
pressure is applied over a small region at the apex to simulate the Im- 
pact of a projectile. Yielding of material In and around the Impact 
zone Is the nonlinear behavior of Interest. Large geometry changes are 
Ignored In the present analysis. 

Figure 2.22 shows the details of the element grids. The square 
panel represents the element grid projection onto the global X-Y plane. 
The middle surface of the shell lies on the surface of a sphere having a 
radius of 60 in. (152.4 cm). The rise along diagonal A-D Is 3.15 In. (8 
cm). Edges C-D and B-D lie on symmetry planes; edges A-C and A-B are 
completely fixed (no translation). A simple von Mises yield criterion 
with associated flow rule is adopted for the nonlinear material 
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Figure 2.22 — Substructures for Non I i 




D 


Top View - Linear Substructure 


Nonlinear Zone 


r Impact Example 



83 


response. The uniaxial stress- strain curve is Idealized as elastic- 
perfectly plastic for simplicity. 

Three types of 3-D solid isoparametric elements are used to model 
the shell. Each node of these elements has 3 translational degrees of 
freedom (u,v,w). In the impact zone, 20 node parabolic elements are em- 
ployed as shown in Fig. 2.22. Two elements are used in the thickness 
direction very near the loaded region. Away from the impact point, one 
16 node thick shell element is adequate through the thickness. A ring 
of 18 node transition elements connect the 20 node and 16 node elements 
without loss of displacement compatibility. A 2x2x2 Gauss Integration 
order is used in all elements. 

Figure 2.22 defines the substructured model components. The outer 
linear region has 483 nodes (1449 DOF), 59 thick shell elements, and 8 
transition elements. Condensation reduces the number of nodes to the 43 
(129 DOF) that Interface with the nonlinear region. The natural node 
numbering scheme places retained nodes last tn DO, which eliminates 
reordering. The highest level structure consists of 56 nonlinear 
parauolic elements In addition to the condensed linear substructure. 
The final nonlinear model has 376 nodes (1128 DOF). The standard model 
for this structure has exactly the same element grid but without sub- 
structuring. This model has 818 nodes (2454 DOF) and 123 elements on 
which the nonlinear solution is performed. The substructured model thus 
reduces the number of DOF by 54$. 
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A total of four load increments were applied to each model with the 
modified Newton-Raphson procedure to distribute residual forces. Itera- 
tions at constant external load were conducted until the Euclidean norm 
of the residual force vector fell below 1$ of the applied load vector 
norm. Tangent stiffness updates were performed before Iterations 2,5, 
and 8 tn each step. Table 2.2 provides the number of yielded Integra- 
tion points, the number of stiffness updates, and the number of Itera- 
tions for each load increment. The fourth load Increment propagated the 
plastic zone Into the third element band from the apex. Sufficient CPU 
timing data was collected during the first four load steps for the 
desired comparisons. 

Figure 2.23 compares CPU time for the standard and substructured 
models. The normalized CPU time used through each load step is ex- 
pressed as a percentage of the total time required for the standard 
model analyzed through load step four. Solution times for load step one 
are nearly equal. The standard model requires slightly more time due to 
the larger number of elements for which strains and stresses are coift- 
puted. The results for step one clearly demonstrate the negligible 
overhead for controlling the substructured solution. The major CPU 
savings with the substructured model begin to occur in load steps three 
and four. At the end of step four, the substructured model used 72$ of 
the CPU required for the standard model. Both solutions appear to have 
reached a steady-state condition (constant slopes) at load step three. 
Linear extrapolation of these curves to one hundred load steps shows 
that solution of the substructured model would require 64$ of the stan- 
dard model solution time. The Improvement from 72$ to 64$ arises from 




Table 2.2 — Computation Summary for Nonlinear Example (Identical Results 


for Substructured and Standard Models) 
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the decreasing significance of equal times in step one and equation 
solver savings In the substructured model. 

The bar graph of CPU distribution among processing modules shown In 
Fig. 2.24 Indicates the source of CPU savings. Input translation and 
strain-stress computation (which includes material updating) consume an 
Insignificant percentage of the CPU time. The equation solver is most 
dominant followed by the stiffness assembler (which Includes element 
computations) and the residual loads generator. Savings in the equation 
solver for the substructured model are approximately 25$ at the end of 
step four. This figure would Increase to 40$ with additional stiffness 
updates during subsequent load steps. The substructured nonlinear re- 
gion with 1128 DOF requires 40$ less CPU to triangulate than the full 
2454 DOF nonlinear model. 

Little savings accrue in the assembler for the substructured model; 
however, this is to be expected. After step one, there is very little 
difference in the operations performed by the stiffness assembler for 
each model. The same number of nonlinear element stiffnesses are 
generated regardless of whether or not the model Is substructured. The 
small savings result from the reduced number of element stiffnesses as- 
sembled In the substructured model. Moreover, this Is an I/O rather 
than CPU Intensive activity. 

Figure 2.25 compares I/O activity for the two models. The I/O ac- 
tivity through each load step Is again normalized by the total I/O 
througn load step four of the standard model. The curves follow the 
same trend as the CPU time In Fig. 2.23. The substructured model 
savings are 33$ at the end of step four, which is slightly better than 
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Figure 2.24 -- CPU Time Distribution for the Example Nonlinear Analysis 




Figure 2.25 — I/O Comparison for the Example Nonlinear Analysis 
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Figure 2.26 — CPU Time Distribution for the Substructured Nonlinear Example 
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the CPU savings. Most of the I/O savings occurred in the equation 
solving file. 

Figure 2.26 provides a breakdown of the CPU time distribution among 
POLO-FINITE subsystems for the substructured solution. A similar com- 
parision for the linear example is given In Fig. 2.20. These two 
figures illustrate the increased percentage of job CPU time expended in 
data management activities for the nonlinear analysis. While it is not 
strictly valid to compare a linear analysis of one structure with the 
nonlinear analysis of another structure, the same trend of Increased 
data management time has been observed In other analyses performed with 
POLO-FINITE. For most linear analyses, data management activities 
represent 10-2556 of the total CPU time. In contrast, data management 
activities consistently require a larger percentage of the CPU time for 
nonlinear analyses. The increase Is attributed to the much larger data 
base sizes (generally a factor of 10), the larger number of data struc- 
tures that must be accessed to obtain nonlinear data, and the large in- 
crease in solution logic that requires movement and duplication of data 
without computation. 
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CHAPTER 3 


DYNAMIC REDUCTION OF STIFFNESS AND MASS MATRICES 


3.1 General 

As demonstrated In Chapter 2, multilevel substructuring provides an 
economical approach to static analysis of very large linear and non- 
linear structural models. Size reduction of the substructures, through 
nodal condensation, yields an exact and economical solution to the 
statically loaded problem. Since dynamic analysis of a finite element 
model requires significantly more computational effort than a static 
analysis of the same model, an analogous reduction scheme would be use- 
ful In dynamics. 

As an Illustration of the usefulness of substructured models in 
dynamics, consider the following example. Assume that a particular 
structural model contains 5000 DOF and has a half-bandwidth of 500. 
Computation of the 50 lowest natural frequencies and corresponding mode 
shapes by a method suitable to the problem characteristics requires 
roughly 7.6(10 9 ) operations. Now suppose that the model can be divided 
Into five Identical substructures, each containing 1200 DOF. Reduction 
of a substructure to 100 Independent DOF while retaining the 10 lowest 

O 

natural frequencies and mode shapes requires roughly 7.5(10 ) opera- 
tions. Since all five of the substructures are Identical, the reduction 
must be performed only once. Assembly of the substructures into final 
form results In the reduced model containing only 300 Independent DOF. 
Since the equations for the substructured model are fully populated, a 
different procedure may be appropriate for computation of natural f re- 
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quencies. Computation of the 50 lowest frequencies and mode shapes 
would require only 2-3(1 0 7 ) operations. Thus, a savings in required 
operations of a full order of magnitude can be realized by using sub- 
structuring on this hypothetical model. These operation counts are, of 
course, highly dependent upon the algorithms used and the model being 
■analysed, but the computational savings is well illustrated. Additional 
savings can be gained in the solution of the equations of motion and in 
the recovery of substructure displacements, strains, and stresses. 

The goal of dynamic reduction methods Is to generate stiffness and 
mass matrices that accurately represent the stiffness and inertia 

characteristics of the substructure with the minimum number of DOh. As 
previously stated, reduction in static analysis Is exact and can be 
mathematical ly viewed as an equation solving technique. In dynamic 
analysis, however, exact dynamic reduction of an individual substructure 
Is dependent upon the unknown frequencies of the total structural 
system. Since these system frequencies are actually objectives of the 
analysis, the analyst must use reduction methods which are either itera- 
tive or frequency independent (and therefore approximate). 

Two classes of methods for dynamic reduction have evolved for use 
with the FEM. The first class, known as Guyan reduction [3.6j, Is an 
extension of static condensation. It Is currently the technique most 
widely used to reduce the number of DOF prior to frequency or transient 
analysis of standard (nonsubstructured) finite element models. The 
method involves elimination of DOF that are assumed to have a negligible 
effect on mode shapes and thus vibration response of the structure. 
Dynamic results, especially strains and stresses, are generally quite 
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sensitive to the choice of DOF to be eliminated. Although the method 
has proven useful for smaller models, its extension to multilevel sub- 
structuring is expected to be only marginally successful. 

The second class of dynamic reduction techniques, termed modal syn- 
thesis, contains methods that rely on a Ray leigh-Ritz transformation of 
■each substructure^ geometric coordinates to a smaller set of 
generalized coordinates. This transformation is usually based on the 
natural frequencies and mode shapes of the Isolated substructures. The 
reduced stiffness and mass matrices for the substructure are defined in 
terms of the generalized coordinates. Assembly of the stiffness and 
mass matrices for the next higher level structure is based on displace- 
ment compatibility. The reduction process can then be successively 
repeated for each additional level of substructur ing. Many techniques 
in the category of modal synthesis have been devised in an attempt to 
select the best combination of substructure modes and displacement com- 
patibility conditions. 

This chapter provides detailed descriptions of the initial formula- 
tions of Guyan reduction and modal synthesis. Derivations of the 
governing equations of the methods are followed by a brief review of 
their respective extensions and modifications. Although many more of 

these alternate techniques have been proposed In the open literature 
than are presented here, those discussed effectively encompass the 
breadth of the topic. The chapter concludes with an evaluation and com- 
parison of the dynamic reduction techniques and recommendations for Im- 
plementation in a general purpose software system. 
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Consider an isolated substructure consisting of simple elements, 
such as the FIN substructure used In the exhaust duct example problem of 
Chapter 2 (see Figure 3.1). Let the Internal boundary of the substruc- 
ture Identify Its interface with other substructures while an external 
boundary corresponds to the physical boundary of the entire structure. 
The undamped, free vibration equation of the substructure, partitioned 
to seperate master (m) and slave (s) DOF, is 



(3.1) 


Master DOF are those that will remain after condensation and are usually 
chosen to lie on the Internal boundary of the substructure. They are 
used for connectivity to adjacent substructures. The slave DOF are 
those to be eliminated and usually lie In the Interior of the substruc- 
ture or on Its external boundary. The natural frequency oj. Is that of 
the complete structural system, not just the Isolated substructure. The 
presence of nonzero off-diagonal blocks CM sm Il and [M nns ]| In Eq. (3.1) Im- 
plies the use of a consistent mass formulation. When a lumped mass 
model Is used, the mass matrix Is diagonal. 

The lower half of Eq. (3.1) can be expanded to 

([K sm ] - u)f[M sm ]) {u m } + (CK SS ] - co?CM SS 3) (u s ) = {0}. (3.2) 

Solving for {u s } In terms of {u m } yields a coordinate transformation 




Internal Boundary (Contains Master C 
External Boundary (Contains Slave DC 


Internal and External Boundaries for Substructure FIN 
from Chapter 2 Exhaust Duct Example 
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which Is dependent on the unknown system vibration frequency co.. If It 
Is assumed that the Inertia forces on the slave DOF are small compared 
to the static forces, the former can be neglected. Thus, the frequency 
dependence Is eliminated and Eq. (3.2) simplifies to 

[K ss ] {u s } - -[K 5m 2 {u m }. (3.3) 

IT) S 

Defining the coordinate transformation CTqD from (u } to (u } as 

{u S } - CTq] tu m } (3.4) 

{u s } can be eliminated from Eq. (3.3) resulting in 

[K ss ] [T^ = -CK sm ]. (3.5) 

As discussed In Chapter 2, [Tq] is evaluated by standard equation 
solving techniques. It Is important to recall that the columns of the 
transformation matrix are composed of the previously defined static con- 
straint modes. The complete substructure displacement vector can be re- 
expressed, using Eq. (3.4), as 
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The substructure strain and kinetic energies are given by 
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(3.7a) 


(3.7b) 


where (u > Is the first time derivative of {u}. These expressions are 
rewritten in terms of the master DOF by substituting Eq. (3.6) Into 
Eq. (3.7) o The resulting Guyan reduced stiffness [K^] and mass are 
given by 

CKj - W™] + [k" s 3 [T s 3 and (3.8) 

- DTh + Ct g 3 t [m ss ] HTg] + Ct/em 5 "’] + Di ms lT^. (3.9) 

For the simpler case of a lumped mass model. Eq. (3.9) reduces to 

cm g ] = Df m ] + CT G a T i;M ss n ct g 3. (3.io) 


In Guyan reduction, the Inertia effects of the slave DOF are not lost, 
instead, the contribution to the mass of these eliminated DOF is 
distributed to the master DOF. It Is assumed that the dynamic response 
of the slave DOF of the substructure Is adequately approximated from 
that of the master DOF by linear combinations of the static constraint 
modes. Note that regardless of which mass matrix formulation is used, 
consistent or lumped, the reduced mass matrix, [Mg]* will be fully 
populated. This situation must be considered when choosing algorithms 
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for solving the free vibration and transient response problems for the 
assembled structure. 

When loads are applied to the slave DOF, they too must be con- 
densed. If the substructure Is subjected to an arbitrary virtual 
displacement, {6u>, the work done by the substructure forces {P} is 

<5W = {6uf {P}. (3.11) 

The statically equivalent condensed forces, (F), applied at the master 
DOF must do the same work during a virtual displacement consistent with 
{6u }, so 

{Su m } T {F} = {<$u} T {P}. (3.12) 

Recalling Eq. (3.6), the condensed force vector becomes 

T 

{P}. (3.13) 

Each substructure has Its stiffness, mass, and loads similarly par- 
titioned and reduced. Assembly of both the reduced substructure mass 
and stiffness into the next higher level follows the procedures outlined 
in Chapter 2. Geometric compatibility between substructures Is 
automatically assured by the use of the master DOF as generalized 
coord l nates. 

The extension of Guyan reduction to multilevel substructuring is 
straightforward. Referring to the terminology of Chapter 2, assume that 
all substructures at level "I" have been assembled either from simple 
elements or level "1+1" substructures. The level W i-1 M substructures 




are built up by selecting slave and master DOF for each substructure at 
level "I", condensing these substructures using Eq. (3.8) and Eq. (3.9), 
and assembling as described above. When the highest level structure is 
reached, the final condensed stiffness and mass matrices can be used to 
form the equations of motion for the entire structure. 

After a free vibration problem has been solved for the highest 
level structure, it may be necessary to recover the portion of the 
system mode shapes contained within lower level condensed substructures. 
This is achieved by simply applying Eq. (3.4) recursively to each sub- 
structure to recover the components of the vibration modes at the slave 
DOF from the master DOF. Recovery of the displacement patterns after a 
transient analysis Is not as elementary. A procedure analogous to the 
computation of partial slave displacements discussed In Chapter 2 must 
be adopted. During the condensation process, loads applied at the slave 
DOF are transformed Into work-equivalent forces at the master DOF. 
Thus, application of Eq. (3.4) to displacement patterns at the master 
DOF does not yield the total response of the slave DOF. A separate 
transient analysis of the substructure with Its master DOF fixed is re- 
quired. The displacements of the slave DOF from this analysis must then 
be superposed with those from Eq. (3.4) to obtain the total response of 
the substructure. Of course, this procedure Is unnecessary when the In- 
ternal response of condensed substructures is of no interest in the 
analysis of the structure. 
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3 .2 .2 Automatic Se I ecti on ai Master IXLE 

To Insure complete geometric compatibility, all DOF at the Internal 
boundaries of a substructure must be included In the set of master DOF. 
However, the set need not be limited to the DOF on the Internal boundary 
of the substructure. Other substructure DOF are possible candidates for 
retention. 

Henshel I and Ong [3.73 have suggested a simple method for 
automating the process of selecting master DOF. The method Is based on 
an assumption fundamental to the development of Guyan reduction; that 
the mass terms corresponding to the slave DOF have a negligible effect 
on the mode shapes. This can be rephrased by saying that *< ss /m ss fs 
large relative to k /m for any pair of slave and master DOF. Thus 
the obvious candidates for master DOF, In addition to those on the sub- 
structure Internal boundary, are those with the smallest ratios k../m. 
1=1, 2,..., n; n=number of substructure Internal DOF. 

The method presented In Ref. [3.73 was oriented towards standard 
finite element models. Its extension to multilevel substructuring re- 
quires only the addition of the internal boundary DOF In the set of 
master DOF. 

The above decision criterion Is equivalent to preserving the lowest 
vibration modes of the substructure. The ratio k../m.. is Interpreted 

i i i i 

as the square of the vibration frequency associated with the "I” th DOF 
while all other substructure DOF are held fixed. Retention of the DOF 
with the smallest values of k../m.. has shown to produce better accuracy 
In the longer modes than random selection of master DOF. 
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3.2.3 Improved Displ acement Recovery 

The frequency dependent transformation from {u m > to {u s }, when cast 
In a modified form. Is useful for Improving recovery of substructure 
mode shapes and displacements. The technique was initially presented by 
Kidder [3.153 and was later reintroduced with the addition of some 
numerical results by Miller [3.203. The technique Is applied after the 
structure equations have been reduced by standard Guyan reduction. The 
frequency dependent transformation matrix resulting directly from 
Eq. (3.2) is 

[TJ - -( ~<*) 2 [M SS 3 + [K SS 3 f 1 ( co 2 [M Sm 3 - [K Sm 3 ). (3.14) 

Expansion of the Inverse term gives 

(-m 2 cm ss 3 + ck 55 ])" 1 - CK ss r' + - 2 ck s, :t 1 [:>i“ik“ 3" 1 + ... <3.i5> 

Ignoring as small the terms containing to to powers greater than two, the 
transformation becomes 

[T 3 - [K ss 3 _1 (-[K sm 3 + (0 2 ([M sm 3 - [M SS 3[K SS 3 _1 [K Sm 3) ) . (3.16) 

to 

As with [T ^3, [T^3 can be evaluated by equation solving rather than by 
computing [K ss 3""^ . 

The transformation [T^3 replaces [T q 3 In recovering the substruc- 
ture mode shapes and displacements during a transient analysis of the 
uncoupled linear equations of motion. Upon solution of the modal equa- 
tions at the highest level, mode shape recovery is achieved via 
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{u s } = [T ]{u m }. (3.17) 

0 ) 

This transformation Is performed for each modal frequency considered so 
as to yield Individual mode shapes within the substructure. Displace- 
ment recovery follows the procedure outlined In Sec. 3.2.1 with the 
total displacement vector equal to the sum of Its modal components. 

The Improved displacement recovery technique has no effect on the 
computed system frequencies and highest level structure mode shapes. 
The argument for Its use Is that frequencies determined from substruc- 
ture mass and stiffness matrices computed with CTqI! are generally 
realistic. However, Improvements are needed In the substructure mode 
shapes from which strains and stresses are derived. 

Miller’s numerical results demonstrated greatly Improved mode shape 
vectors over standard Guyan reduction for frame structures. However, 
the effectiveness of the method Is limited to the lower frequencies of 
vibration. This Is because the size of the truncated terms of the 
series In Eq. (3.15), and thus the error In CT ], grows with increasing 
values of u>. 

In order to use Improved displacement recovery, the equations of 
motion of the highest level structure must be solved In their uncoupled 
form so that modal vibration frequencies can be computed. When a tran- 
sient analysis of coupled equations Is performed, as In nonlinear 
analysis, the displacements are computed directly and thus, this method 
Is not appl I cable. 
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3.2.4 Eva I uatlon of Guyan Reduction Techniques 

Guyan reduction techniques have yet to be applied to multilevel 
substructured models. Therefore, their performance in such an applica- 
tion remains unknown. Nevertheless, Guyan reduction has some important 
advantages that could make it an attractive approach to dynamic reduc- 
tion In particular cases. First, its development Is taken directly from 
the static condensation approach. This allows the reduction scheme to 
be readily added to existing software which Is currently capable of 
handling multilevel substructured models for static analysis. Secondly, 
reasonably good numerical results have been achieved with Guyan reduc- 
tion In computing system frequencies for small models. This feature 
makes the technique attractive for preliminary vibration analysis. 
Lastly, Guyan reduction is the least expensive of al I dynamic reduction 
techniques. 

Some drawbacks of Guyan reduction are evident and must be con- 
sidered. The success of the method Is highly dependent upon the choice 
of master DOF. This Is further complicated by the need to include in- 
ternal boundary DOF in the set of master DOF for substructured models. 
The result is likely to be a decrease in the degree of reduction capable 
for large, multilevel substructured models. Another problem Is the 
quality of the mode shapes. Accurate prediction of strains and stresses 
requires that displacement vectors and mode shapes be well formed. The 
ability to achieve this goal with Guyan reduction Is still in doubt. 
However, there are circumstances in which it is not necessary to recover 
strains and stresses within a condensed substructure. In such cases, 
i.e. when the response of only the highest level structure Is of major 
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concern, Guyan reduction may perform quite well. 

Guyan reduction techniques provide an important first step in the 
development of more sophisticated dynamic reduction methods. As will be 
seen later, the procedures derived above actually represent a 
degenerated case of modal synthesis. Guyan reduction does not appear 
adequate for general application to dynamic analysis of multilevel sub- 
structured models. However, with automatic selection of master DOF and 
improved displacement recovery, It does hold potential as an economical 
approach to dynamic reduction of certain models. 
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3.3 Moda 1 Synthesis 
3.3.1 I ntroduction 

Modal synthesis was developed expressly for use with substructured 
models. Although this discussion Is limited to reduction of finite ele- 
ment models, the techniques have also been applied to distributed 
systems. Ail modal synthesis transformations are based on Ray leigh-Ritz 
arguments. The procedure Involves the derivation of a transformation 
matrix composed of a truncated set of mode shape vectors that adequately 
describe the dynamic characteristics of the substructure. This set is 
fewer in number than the number of independent DOF contained In the sub- 
structure. With the transformation matrix, individual substructure mass 
and stiffness matrices are converted from geometric coordinates Into a 
reduced set of generalized coordinates. The generalized mass and stiff- 
ness matrices for each substructure are then synthesized while main- 
taining geometric compatibility along internal boundaries to form 
similar matrices for the next higher level structure. In a multilevel 
substructured environment the transformation and assembly processes are 
performed recursively at each level. 

There are two basic operations that must be performed with any ap- 
plication of modal synthesis. First, an approach must be chosen for 
selecting the set of substructure mode shapes from which the reduced 
substructure matrices are computed. Second, a procedure Is needed to 
enforce geometric compatibility along substructure internal boundaries. 
The numerous modal synthesis techniques proposed in the literature vary 
in how these concepts are Implemented. 
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The Initial formulation of modal synthesis, accredited to Hurty 
C3.11, 3.12], has been extensively modified and enhanced. The following 
section presents Hurty’s formulation, commonly referred to as the fixed- 
interface method, while the later sections describe suggested improve- 
ments or alternatives to the method. 

3.3.2 Ft xed- Interface M albfld. 

The origin of modal synthesis techniques lies in Hurty's fixed- 
interface method. However, details of the procedure presented here 
parallel the development by Craig and Bampton [3.3]. Hurty's develop- 
ment required a distinction between "statically determinate" and "redun- 
dant" constraints. Statically determinate and redundant constraints are 
artificial constraints applied at all master DOF. They are imposed in- 
dependently of any actual physical constraints that may already exist on 
the substructure at slave DOF. The set of statically determinate con- 
straints serve to restrain any rigtd-body motion that may be possible. 
Redundant constraints are those applied to the remaining master DOF. 
Figure 3.2 illustrates the application of these constraints to a simple 
two-dimensional plate. Three statically determinate constraints are 
needed to restrain translation and rotation. The five remaining master 
DOF have redundant constraints applied. Craig and Bampton treat the two 
sets of constraints simply as boundary constraints. 

Consider again the undamped, free vibration equation for an 
isolated substructure composed of simple elements and partitioned to 
separate master and slave DOF 



i 09 



• Master Node (2 DOF/NODE) 
o Slave Node (2 DOF/NODE) 

Statically Determinate Constraint 
JlL Redundant Constraint 

SJ77777 


Figure 3.2 — Statically Determinate and Redundant Constraints on a 


Simple Two Dimensional Substructure 
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As In Guyan reduction, a static transformation from the master to the 
slave DOF can be written 


{u S } = U c ]{u m }. (3.19) 

As with simple Guyan reduction, the set of master DOF In modal synthesis 
methods Is not limited to those DOF on Internal boundaries. DOF In the 
Interior of the substructure may be retained as well and may possibly 
Improve the solution. Inspection of Eq. (3.2-3. 4) reveals that Dj> c I] Is 
Identical to [T^, the static constraint modes. For consistency In the 
following development, [cf) C I] Is used to represent the static constraint 
modes. 

If the set of master DOF Is restrained from displacement, 
Eq. (3.18) reduces to 

CK SS ]{U S } - Wj[M ss ]{u s } = {0}. (3.20) 


The solution of this eigenvalue problem yields the matrix of fixed-fixed 
normal modes of vibration, having the same order as [K ss ] and 
CM bS 3. The computed vibration frequencies, w. , are those of the 
Isolated substructure. 

To reduce the substructure mass and stiffness matrices, a transfor- 



ma+ion to generalized coordinates, {q }, Is defined as 



(3.21) 


The f I xed- i nterf ace transformation, [T^J, Is derived from the static 
constraint and the normal modes as 


Ct f : 



(3.22) 


in which C<j> n U is a rectangular matrix of retained modes from C(j> n []. In. 
general, the modes corresponding to the lowest natural frequencies, w , 
are retained in [q n H. The slave displacements, {u s }, are now dependent 
on both the static constraint modes and the retained normal modes of the 
isolated substructure. 

Two observations regarding Eq. (3.21) are noteworthy. First, the 
generalized coordinate subvector, {q m }, corresponds precisely to the 
master set of geometric coordinates, {u m }. This proves useful in 
guaranteeing geometric compatibility with adjacent substructures. 
Secondly, as the number of mode shapes In C$ n 3 is reduced, the transfor- 
mation shrinks to just the static constraint modes and thus, the fixed- 
interface method degenerates to Guyan reduction. 

The reduced stiffness and mass matrices in generalized coordinates 
are obtained by maintaining equivalence of kinetic and potential ener- 
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gies between the two coordinate systems. The resulting forms are 


ck f : = ctjtkxt.: 
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[m f : = [t f ]'cmxt f : = 


:m g d 


cr]V s x^ c ] 

+ i:$ n ] T CM sm : 


[M mS ]K n ] T + 
^ c : T CM ss x^ n : 
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(3.24) 


For a lumped mass formulation of the substructure 


[m f :i 


cm g ] 


z* c i T o? s xri 


[^ n : T [M ss x<i) c ] ! i n 


(3.25) 


[K H and are the Guyan reduced stiffness and mass matrices. The 

G G 

-2 

identity submatrix in [M f 3 and the submatrix Cw H in CKp-H result from 
the orthonormality of C$ n 3. is a diagonal matrix of natural fre- 

quencies corresponding to the modes retained in [<j> n 3. 

Reduction of substructure loads, {P}, follows the same virtual work 
argument used In Guyan reduction. The resulting generalized force vec- 
tor, {F}, is 

{F} = CT f ] T {P}. 


(3.26) 



Before proceeding, a simplified notation is introduced. Let 


[K p ] = 

with the relationships to Eq. (3.23) and Eq. (3.24) established by in- 
spection. 

Although assembly of the reduced substructure stiffness and mass 
matrices Is routine, an illustration of their form is useful. For an 
assembly of "r” substructures with the master DOF entered last 


(3.28a) 
(3.28b) 

Since the master DOF do not participate in the normal modes, no coupling 
between substructures exists outside the submatrices CK mm ] and [M mni ], 
which are the assembled Guyan stiffness and mass. From a data storage 
and computational viewpoint, these forms for DO and CM] are ideally 
suited to hypermatrix methods. 

The synthesis process for one level of substructur i ng is now com- 

¥r 

plete. if [K] and [MU are the stiffness and mass matrices for the 
highest level structure, the differential equation of motion can be 
written and solution for displacements can proceed. If the highest 
level structure has not yet been reached, [K] and QG are reduced just 
like any other substructure. 
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In summary, the f ixed- interface method employs static constraint 
modes and a truncated set of fixed-fixed normal modes to achieve a 
reduction In the number of independent substructure DOF. Geometric 

coordinates at internal boundaries are retained in the set of 

general I zed coordinates to assure displacement compatibility between 
substructures during assembly. This mixture of geometric and 
generalized coordinates In the substructure equations has resulted In no 
reports of numerical difficulties and none are expected. 

3.3.3 Free- Interface Method with Interface Loadl ng 

In the free- interface method of modal synthesis [3.5, 3.10H, the 

transformation to a reduced set of generalized coordinates relies on a 
truncated set of free-free normal modes. These mode shapes are computed 
for the Isolated substructure with geometric constraints applied only as 
they occur In the actual structure. Internal boundary (master) DOF are 

not arbitrarily Identified and fixed. This approach allows rigid-body 

modes to appear in the set of substructure normal modes unless suf- 
ficient physical constraints actually exist. The transformation to 
generalized coordinates also neglects the static constraint modes common 
to both Guyan reduction and the f Ixed- interface method. The transformed- 
tlon is simply 

{u} == U n ]{q} (3.29) 

where [$ n 3 contains the truncated set of free-free normal modes cor- 
responding to the lowest substructure natural frequencies. 



The development by Goldman £ 3.53 Includes a distinction between 
rigid-body modes and free-free elastic modes which is not required in 
Hou’s method £3.103. These two authors use similar approachs to sub- 
structure assembly. Geometric compatibility between substructures Is 
enforced by writing equations of constraint In modal coordinates for In- 
ternal boundary DOF. These constraint equations are used to combine the 
generalized displacement vectors for all Isolated substructures Into a 
single generalized displacement vector for the final structure. The 
corresponding reduced stiffness and mass matrices are thus generated. 

A significant difficulty arises with the extension of the free- 
interface method to multilevel substructuring. In achieving one level 
of substructuring, all geometric coordinates are transformed to 
generalized coordinates. Further substructuring is complicated by the 
absence of geometric coordinates which are useful In assuring continuity 
of displacements at substructure Internal boundaries. Beyond the lowest 
level. It will be necessary to modify the process of developing con- 
straint equations to link together the substructures. Rather than 
defining the constraints In terms of the generalized coordinates at the 
level being assembled, they must be written In terms of the geometric 
coordinates at the lowest level. Such a task has yet to be Investigated 
and the requirements for Its Implementation remain unknown. 

The mode shapes of a substructure with free (or fixed) boundaries 
are not totally representative of the substructure’s response in the as- 
sembled structure. This follows because the stiffness and Inertia ef- 
fects of the adjoining substructures have not been included In computing 
the dynamic modes of the isolated substructure. Interface loading £3.2, 



3.9] is a technique which incorporates some of these effects in an at- 
tempt to make the isolated substructure modes more like the modes for 
the entire structural system when the free- interface method is used. 
The approach is to modify the stiffness and mass of each substructure 
prior to extracting its free- free normal modes. 

Let the DOF of the substructure under consideration reside In set A 
and let all remaining DOF In the structure reside in set B. In general, 
set B may contain more than one substructure. Upon assembly, set A will 
Intersect set B over the subset AB. Consider subset AB as the master 
DOF for both sets A and B. Figure 3.3 illustrates this approach as ap- 
plied to the shallow shell example of Chapter 2. The substructure under 
consideration Is the linear zone and it forms set A. The remainder of 
the structure, the nonlinear zone, falls in set B. Subset AB is the in- 
terface between the two sets. 

By the same transformation as used in Guyan reduction, the 
displacements of DOF in set B can be written in terms of those In subset 
AB. Recognizing that subset AB also defines the master DOF of set A, 
the displacements of DOF in set B are effectively expressed in terms of 
the master DOF of set A. From these relationships, the interface loaded 
stiffness and mass for set A are derived. The effect of the Interface 
loading is to add the Guyan reduced stiffness and mass of set B to the 
Internal boundary terms (master DOF) of A's stiffness and mass. 

The free- free normal modes of the substructure in set A are com- 
puted using the modified stiffness and mass matrices. Each substructure 
that will be reduced is Identified as set A with the remaining substruc- 
tures lumped into set B. A new Interface loading effect is computed and 




Figure 3.3 — DOF Sets for Interface Loading of the Linear Substructure 
from the Chapter 2 Shell Example 




the free- free analysis performed. Actual geometric constraints are ap- 
plied as they occur in either set. After solution for the normal modes 
of each substructure is complete, the substructures are assembled using 
the truncated mode sets and the original, unmodified substructure stiff- 
ness and mass matrices. The interface loaded stiffness and mass are not 
used in the assembly process. Their use Is limited to computing free- 
free normal modes for the substructures. 

3.3.4 Branch Mode Anal y sis 

Branch mode analysis [3.4, 3.9] is a hybrid method of modal syn- 
thesis which Incorporates features of both the fixed- interface and' the 
free- interface methods. The following procedure is one of several that 
fall under the class of branch mode analysis [3.13]. 

In branch mode analysis, one substructure is selected as the main 
body and all adjacient substructures are designated as branches. The 
analysis starts with a determination of the free-free component modes of 
the main body. Interface loading may be employed prior to computing the 
free- free modes. The topology of the main body is then expanded to in- 
clude all substructures on Its boundaries. This is done by writing the 
displacements for the main body and its adjacent branches in terms of 
the free-free modes of the main body and the fixed-fixed normal modes 
from each branch. Using the transformation to modal coordinates In the 
structure kinetic and potential energy equations, the reduced branch 
stiffness and mass matrices are computed. 



When substructures are remote from the initially chosen main body, 
a number of solution schemes may be used. One approach is to Identify 
one main body and its branches and perform the foregoing reduction. 
This reduced substructure group is then treated as a branch to a newly 
selected main body. Two more ei genprob I ems are solved and the process 
Is repeated until each substructure has been joined. 

3.3.5 D ynamic Stillness. Matrix 

An iterative method for obtaining an exact solution to the frequen- 
cy dependent dynamic reduction problem has been developed by Leung 
C3.173. The method Is closely related to the fixed- interface method but 
includes the frequency dependence In the transformation to the reduced 
set of coordinates. The first step of the method is to partition the 
stiffness and mass of each substructure and compute the Guyan reduced 
matrices [K G H and [M^U* With the interna! boundary DOF identified as 

master DOF, the fixed-fixed normal modes, (4> n }» and associated frequen- 

_2 

cfes, uk, are computed for each substructure. The exact substructure 
dynamic stiffness matrix, CD], can then be derived as 

CD] - Ck g ] - a 2 CM G ] - « 2 [G]CAIG] T (3.30) 

2 

where = an unknown system frequency 

CG] - [M mS I<j> n ] - CK ms I* n I5 2 r 1 
-2 

Coj. 3 = a diagonal matrix of substructure frequencies 
1 2 
Cx3 - a diagonal matrix of elements “? ~i 

co - 0 . 

With no effect on the order of [D], all normal modes in Dj> n IS may be used 
or the set may be truncated to Include only those associated with the 
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lowest frequencies. 

Using a theorem stating that the dynamic mass matrix, [M(w)], 

_2 

equals the partial derivative of CDH with respect to oj , the substruc- 
ture reduced mass Is given by 

= CMgU + CGXQXG] 1 (3.31) 

-4 

00 . 

where [Q] = a diagonal matrix of elements 1 ■= ~-~- 

(w - err 

The matrices CK^] and [M(w)3 for each substructure at the same 
level are assembled by the method used In Guyan reduction to obtain syn- 
thesized stiffness and mass matrices for the next higher level. The 
process Is then repeated until the highest level structure is reached. 

The solution process for this formulation requires Iteration at 
each level of substructures. The unknown system frequency, ft, must be 
Initially estimated and then Iteratively Improved until convergence Is 
attained. One Iteration Involves building [M(ft)H and [DU for each sub- 
structure until the highest level structure Is reached. Then the eigen- 
problem at the highest level Is solved, 

= {0}. (3.32) 

From this solution, the next value of ft Is chosen for use In the fol- 
lowing iteration. Although the fixed- fixed normal modes need be com- 
puted only once for each substructure, each Iteration requires that the 
mass matrix for the intermediate level substructures be recondensed. 



It Is stated that a partial vibration case results when Q = oj. for 
a substructure and thus drives the dynamic mass to infinity. However, 
there are no clear guidelines given for handling this potential in- 
stability in the method. 

The method is shown to yield good accuracy in computing system fre- 
quencies but results for mode shapes are not given. Also, no comparison 
of computational efficiency is made with the more standard methods of 
modal synthesis. 

3.3.6 Attachment Modes and Interface Mode Sets 

Bamford, Wada, Garba, and Chisholm C3.il] Introduced the concept of 
attachment modes as an additional set of static modes used In modal syn- 
thesis. An attachment mode deftnes the response of the substructure to 
a unit force applied at an internal boundary DOF while all other inter- 
nal boundary DOF remain free. The motivation for using attachment modes 
is that their use can be expected to reduce the number of normal modes 
necessary to accurately describe the displacement behavior of the sub- 
structure. 

In the transformation from generalized to geometric coordinates for 
the substructure, attachment modes are combined with static constraint 
and normal modes. The normal modes may be computed with Internal boun- 
dary DOF fixed, free, or a combination of the two. Bamford* s transfor- 
mation matrix Includes a set of rigid-body modes, however it is noted 
that no distinction between rigid-body modes and static constraint modes 
is required. A drawback In the use of attachment modes is that they are 
not always linearly independent of the static constraint and normal 



122 


modes. 

In an effort to establish a more rational approach to substructure 
synthesis, HIntz [3.83 grouped combinations of the four mode classes: 
rigid-body, static constraint, normal, and attachment. Into five dif- 
ferent Interface mode sets. The four mode classes are Illustrated in 
Figure 3.4 as they apply to the FIN substructure of Chapter 2. Each of 
the five sets Is claimed to be complete In that It precisely represents 
the original finite element model for static and dynamic response. Im- 
plications of truncating a selected interface mode set are considered 
and thus guidelines are developed for retaining accuracy with a reduced 
size model. Since each of the five mode sets represents the substruc- 
ture differently, truncation of each set has varying Impacts on the sub- 
structure model. The guidelines are directed toward establishing a mode 
set that will allow maximum truncation of the normal modes while 
retaining the detail of the original finite element model with respect 
to statically imposed interface forces and displacements. 

In numerical results presented by HIntz, the fixed- Interface method 
(defined by one of the five Interface mode sets) gave good results over 
a broader range of frequencies than did those methods using attachment 
modes. However, use of attachment modes did produce more accurate 
vibration frequencies In the low frequency range. 
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3.3.7 Improved Displacement Recovery i n Moda I Synthesis 

The frequency dependent transformation previously discussed with 
regard to Guyan reduction, Eq. (3.14), has been applied in a limited 
fashion to modal synthesis by Kuhar and Stahle C3.16H. The method con- 
sists of reducing the order of the eigenvalue problem for the highest 
level structure after modal synthesis has been performed. If ft is a 
frequency about which the eigenproblem for the highest level is to be 
reduced, the synthesized matrices can be reduced with the transformation 

[T] = -([K ss ] - ft 2 [M ss ])~ 1 (0< Sm ] - ft 2 [M sm ]). (3.33) 

The reduced matrices form an eigenvalue problem of smaller size than 
that defined by the synthesized matrices. 

As with Guyan reduction, the eigenproblem is evaluated more ac- 
curately at the frequencies near the reduction value, ft. (In Guyan 
reduction, the reduction frequency Is ft = 0.) Thus, the choice of master 
DOF for the highest level structure should be directed towards 
preserving the frequency content of the model around ft. Although selec- 
tion of the master DOF is Initially a matter of Judgement, It can be Im- 
proved In the second and later Iterations If necessary. If the frequen- 
cy range of Interest is broad, it may be advantageous to re-sol ve the 
problem for different values of ft. 

Upon solution of the reduced eigenproblem at the highest level, the 
set of mode shapes, d4> IT> ] , will be In terms of the master DOF. If a com- 
puted frequency w. Is equal to the reduction frequency ft, then the cor- 



responding mode shape is obtained exactly by 


4J 




Likewise, the full set of mode shapes Is obtained from 



(3.34) 


(3.35) 


but the vectors in C$13 are only approximate. Based on the coefficients 
In the mode shape vectors corresponding to frequencies within the range 
of Interest, new master DOF can be identified to improve the accuracy of 
the reduction. 

3.3.8 Res] dua l Mass aM Easi-dual Flexibility 

A hybrid method of substructure reduction with the capability of 
modeling a more general condition of internal boundary constraint was 
developed by MacNeal [3. 18] and implemented to a limited extent in 
NASTRAN. The normal modes of the substructure are computed with inter- 
nal boundary DOF free, fixed, or a combination of the two. The decision 
regarding boundary DOF fixity is based on the nature of the adjacent 
substructures in the model. If an interface DOF in one substructure is 
held fixed, the corresponding DOF in the adjacent substructure must 
remain free in order to avoid an overconstraint problem at the inter- 
face. This procedure requires that Information on the topology of the 
substructures be provided by the analyst prior to modal analysis of any 
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Individual substructure. MacNeal’s formulation also requires that the 
mass at fixed boundary DOF be distributed to free DOF in the substruc- 
ture since fixed DOF do not participate in the vibration response. 

The representation of the substructure can be improved by including 
static approximations to the effects of the truncated higher modes. The 
improvement is in the form of a residual mass matrix for the synthesis 
Involving fixed-fixed normal modes and a residual flexibility matrix 
when free-free normal modes are used. For the case of hybrid modes 
(some internal boundary DOF fixed, others free), both matrices are for- 
mulated. Since these matrices are obtained by a static derivation 
(w= 0), they are valid only at vibration frequencies that are low com- 
pared to the lowest mode of the substructure. 

Rubin C3.21U proposed an Improvement to MacNeal's residual flex- 
ibility for substructures with free-free normal modes only. The ap- 
proach is shown to have better convergence than the methods of MacNeal 
and Hurty but like MacNeal’s, its application is strictly limited to the 
low frequency range. Rubin’s method is further restricted to the use of 
f ree- i nterf ace normal modes and rigid- body modes. Thus, MacNeal’s 
residual mass Is not included in the development. 

3.3.9 Low- Order Polynomial lLajjS f , Qrm a.t. i .Q,as. 

One of the primary differences among the reduction techniques 
discussed thus far lies in the definition of the transformation to 
generalized coordinates. Each method that has been reviewed uses some 
combination of static constraint and substructure normal modes. 
Meirovttch and Hale [[3,193 have recognized that It is not necessary to 
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consider only substructure modes for use In the transformation. They 
have shown that It Is sufficient for substructures to be represented by 
admissible functions that are from any complete set. One such set con- 
sidered useful for continuous systems Is a set of low-order polynomials. 
As applied to finite element models, the polynomials are defined over 
the domain of the substructure and satisfy the geometric boundary condi- 
tions. The shape vectors are then built by sampling the polynomials at 
the spatial coordinates of each substructure DOF. 

The transformation to generalized coordinates Is 

{u} - [T p ]{q}, (3o36) 

where CT p ] Is the set of shape vectors derived from the low-order 
polynomials. As with the free- Interface method, the generalized 
coordinates, {q}, contain no geometric equivalents for substructure as- 
sembly. Compatibility between substructures Is enforced by a weighted 
residual method using spatial Dirac delta functions for continuous 
models or equivalently, unit vectors for discrete models. For a finite 
element model, the effect is achteved by writing equations of constraint 
matching the displacement at the point shared by two substructures. 
Using the constraint equations, a constraint matrix Is produced which Is 
used to link together the otherwise isolated substructures. From this 
operation Is obtained the synthesized stiffness and mass matrices In 
generalized coordinates. 
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A difficulty lies in development of the constraint matrix. In the 
elementary example problems presented, the matrix was built by inspec- 
tion. This is not an acceptable procedure for general application. As 
observed by the authors, an alternative formulation can be defined which 
relieves the problem of constructing the constraint matrix. The ap- 
proach Is simply to use the same transformation matrix as used In the 
fixed- Interface method but with fixed-fixed normal modes replaced by 
lov/-order polynomial shapes. Thus 



(3.37) 


where C<f> P II is the set of low-order polynomials. For this approach, 
polynomials which have zero values at the Internal boundary nodes must 
be chosen to simulate the fixed-fixed boundary conditions. Geometric 
compatibility is enforced by the presence of the static constraint modes 
[<f> c ] and retained internal boundary nodes in geometric coordinates. 

There are two trade-offs that must be made in using the transforma- 
tion of Eq. (3.37) instead of that in Eq. (3.36). The static constraint 
modes must be evaluated causing some additional computational expense. 
Also, the size of the set of admissible polynomial functions is reduced 
by the requirement that the functions be zero- valued at substructure in- 
ternal boundaries. 

In contrast to individual finite elements, the geometry of sub- 
structures is not of a nature that can be easily classified. Thus, 
selection of the appropriate polynomials is a matter of experience and 
Judgement. Because no automated selection method has yet been devised. 



this approach to modal synthesis Is currently not capable of performing 
dynamic reduction In a general finite element system. However, the 
procedure remains one of the more actively pursued techniques and does 
hold potential. 

,3.3.10 Eml.JiatLon at Moda I Synthesis Techniques 


In the study of modal synthesis techniques, several advantageous 
characteristics of an Ideal method can be identified. A discussion of 
these characteristics as evaluation criteria will be useful In comparing 
the various methods presented. 

1 . Efficiency of the Reduction Method 

The efficiency of a dynamic reduction method is In- 
fluenced by a number of factors. First, the method must 
result In an accurate reduction of the substructure stiffness 
and mass. An efficient method will yield synthesized stiff” 
ness and mass matrices that accurately maintain the dynamic 
characteristics of the substructure with the minimum number of 
DOF. Second, the degree of analyst participation should be 
limited to simply the definition of the model and specifica- 
tion of the solution type. A method should be automatic once 
the solution process begins, hence eliminating the need for 
the analyst to Interpret intermediate results and restart the 
process. This Is not to Imply that the analyst should sur- 
render control of the solution process. Instead, he should be 
relieved of the burdensome task of supervising the com- 
putational process. Third, the synthesis method should be ef- 
ficient In Its use of the computer. Given the problem size, 
algorithms should be chosen that minimize the required com- 
puter resources, including CPU time and I/O. The number of 
arithmetic operations performed should be predictable rather 
than dependent upon an arbitrary test for convergence of an 
iterative process. 

2. Ap^lleaMLlly is geneol Prob-ijans. 

A wide variety of dynamics problems exists for which 
modal synthesis Is needed to achieve an economical and ac- 
curate solution. A synthesis method used In a general purpose 
FEM system should be capable of representing substructures 
over a broad range of geometries and with various types of 
boundary constraint. Also helpful would be the ability to In- 
corporate experimental data, such as mode shapes and natural 
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frequencies, into the substructure model. 

By necessity, finite element analyses of nonlinear struc- 
tures are performed incrementally. As the effects of non- 
linear materials and geometry occur, the mode! must be refor- 
mulated to accurately track the true response of the struc- 
ture. Therefore, dynamic reduction methods must lend them- 
selves to this incremental solution process. 

3. I ndependence of I n d f v I dua I Substructures 

Often times the analysis and design responsibilities of 
the various components of a structure are distributed among 
different organizational groups. This separation of respon- 
sibilities has many advantages and should not be encumbered by 
the synthesis method. Therefore, the method used should con- 
sider each substructure as an isolated entity in evaluating 
Its dynamic response prior to system synthesis. Topology of 
the substructures should not be considered until the equations 
at the next higher level are ready for assembly. 

4. Ease of Reanal ysi s 

The most reliable test for convergence of a dynamic 
reduction method is to re- solve the problem with a more highly 
refined model (more Independent DOF). The addition of more 
DOF to the model can be a relatively simple task, achieved at 
little expense, or It could be as difficult and expensive as a 
complete reanalysis of each substructure. The ideal approach 
allows simply the addition of previously neglected terms to 
improve the accuracy of the reduction. These terms generally 
take the form of truncated substructure normal modes. 


5. 


aM 


Accuracy of results Is Important in two respects. Well 
defined modal response data Is needed to accurately synthesize 
the higher level structures for frequency and transient 
analysis. Also, in returning to the lower levels for recovery 
of strains and stresses, the quality of the displacement vec- 
tors is critical. Accurate stresses require that displacement 
gradients be well formed. Intimately tied to accuracy of 
results is the numerical precision with which computations 
must be performed. Operations such as triangul at Ion in equa- 
tion solving can have a significant impact on final accuracy 
and the need for such operations should be considered in 
selecting the reduction method. 

The potential for problems with the stability of opera- 
tions in the reduction methods can often be identified by 
close examination of the development of the methods. Typical 
problem areas are linear dependence of the vectors comprising 
a transformation matrix and the divide- by- zero singularity. 
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With the criteria established, each synthesis method can be 
evaluated. The goal of this evaluation is to isolate one or two methods 
of modal synthesis that will be most useful In the general purpose FEM 
system. 

The f ixed- interface method successfully satisfies four of the a.bove 
five criteria. The method is simple to apply and results In a sig- 
nificant size reduction of properly substructured models. It does not 
require consideration of substructure topology prior to assembly thus 
preserving substructure independence. Repetitive tests for convergence 
are straightforward, requiring simply the addition of more normal modes 
In selected substructures. The entire problem need not be re-solved to 
Incorporate the additional modes. Thus the cost of convergence tests is 
low and predictable. Orthogonality of the transformation matrices en- 
sures stability of the method and accuracy has proven favorable for many 
problems. One drawback of the method Is its limited capability to use 
experimental data. Since all master DOF are fixed during the computa- 
tion of normal modes, the structural components must be tested In a like 
fashion. 

Interface loading has proven necessary to the use of the free- 
interface method. The approximate Inertial effects of adjacent sub- 
structures serve to improve the displacement gradients at the internal 
boundaries, thus yielding more accurate strains and stresses [ 3 . 143 . 
Although interface loading Is helpful in computing accurate stresses, 
its use eliminates the ability to maintain substructure Independence. 
As previously mentioned, extension of the free- interface method to mul- 
tilevel substructuring might be a difficult task to accomplish. 
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Branch mode analysis was originally developed for chain-type struc- 
tures, e.g. piping networks. As such, its application to analysis of 
more general structures is cumbersome. The method requires solution of 
twice as many eigenproblems as the fixed- interface method and it is in- 
timately dependent upon substructure topology. Convergence testing 
using a more refined model becomes impractical as complete reanalysis is 
necessary. 

Leung’s iterative method appears limited to eigenvalue extraction 
over a narrow band of frequencies or analysis of transient response to 
harmonic forcing functions. This is because the system equations must 
be synthesized for each modal frequency considered. For multilevel sub- 
structured models, the need for repetitive condensations is expected to 
make the method prohibitively expensive when the frequency range of In- 
terest is broad. As mentioned In the development of the method, some 
special precautions would also be necessary to maintain stability when 
the system frequency is very close or equal to a substructure modal fre- 
quency. 

Although attachment modes are useful in reducing the number of 
retained substructure normal modes, their potential failings are not 
well illustrated in the literature. It Is noted that attachment modes 
may cause Ill-conditioning In the transformation to generalized 
coordinates C3.1D. What must be recognized is the difficulty with which 
linearly independent attachment modes are chosen. The selection of at- 
tachment DOF for a substructure Is by no means an intuitive process. It 
is expected that even the experienced analyst will require a trial-and- 
error approach to their selection. 
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Efficient methods for extracting a limited number of eigenpairs in 
large systems effectively make the highest level structure reduction 
proposed by Kuhar and Stahle unnecessary. Their process yields a 
reduced size equation of motion via an iterative procedure. The same 
results can easily be achieved by solving for selected modes of the un- 
reduced synthesized equations and then formulating a transformation with 
these modes and appropriate static modes. This transformation can then 
be used to reduce the size of the synthesized equations. 

MacNeal’s residua! mass and flexibility method provides an 
automatic approach to modal synthesis which Includes general boundary 
constraint capabilities. However, the method does have some limita- 
tions. It is correctly argued that fixed coordinates do not participate 
in the vibration behavior of the substructure. What is not considered 
is that in a multilevel substructured model, the coordinates held fixed 
at level ”1+1 ” may be considered free at level ”1". If the mass compo- 
nents of the constrained DOF are redistri buted as suggested, the inertia 
characteristics of those DOF will not be accurately preserved. Also, 
the need to define substructure topology and properties prior to as- 
signing Internal boundary constraints hampers the effort to retain In- 
dependence in computing individual substructure response. Another draw- 
back is the potential ill-conditioning of the residual flexibility 
matrix when many modes are retained in the transformation. The advan- 
tages of using residual mass and flexibility matrices appear marginal in 
that they Improve only the low frequency response and that residua! mass 
is implicitly contained in Hurty f s method [3.18], 
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Polynomial transformation techniques hold tremendous potential as 
dynamic reduction methods. A great savings in computational effort is 
possible by evaluating low-order polynomials rather than computing sub- 
structure vibration modes. Unfortunately though, the approach to 
selecting the polynomials has not yet been standardized. As such, the 
techniques are not Immediately useful for general application In FEM 
software. 

3.4 Select.! on of Methods for Dynamic Reduction 

In view of the above discussion, selection of dynamic reduction 
schemes becomes simple. The two broad types of structural models, those 
with and without substructur I ng, exemplify the need for two distinct 
methods for dynamic reduction. The two methods considered most suitable 
are Guyan reduction and the fixed- Interface method of modal synthesis. 
Perhaps the greatest advantage of both of these methods Is their concep- 
tual simplicity. This point cannot be overemphasized. It Is Imperative 
that the analyst have a complete understanding of the analysis methods 
he uses. Without this understanding the chances of achieving a meaning- 
ful dynamic analysis are remote. The obvious disadvantage of this 
philosophy Is the sacrifice In sophistication that the experienced 
analyst makes if he is to use the general system. One might expect, 
however, that special purpose FEM systems could be used when the need 
for more refined techniques is not satisfied elsewhere. 
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CHAPTER 4 


COMPUTATIONAL ALGORITHMS FOR DYNAMIC ANALYSIS 

4.1 General 

The solution of structural dynamics problems by the FEM requires 
computational capabilities which are not necessary In static analysis. 
The two most Important of these. In terms of computational efficiency, 
are elgenproblem solution and transient response analysis. These two 
operations form the core of the dynamic analysis process. Their proper 
Implementation and use is essential to the success of the dynamic 
analysis. Other less Important features include mass matrix and damping 
matrix formulation and the use of experimental data. Little com- 
putational effort is expended on these operations but their use adds 
generality and completeness to the solution strategy. 

Elgenproblem solution is the single most expensive operation in the 
modal synthesis process. While the approach to modal synthesis will 
control the quality of the solution, eigenvalue analysis can be expected 
to control the cost. An understanding of the various methods for eigen- 
value analysis is, therefore, necessary if the analyst is to achieve an 
economical and accurate solution. 

Solution of the differential equations of motion yields the com- 
plete structural response to the transient loading. This solution can 
be obtained by either of two approaches: mode superposition for linear 
systems and time-history Integration for both linear and nonlinear 
systems. 
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The remainder of this chapter is divided into three sections. The 
first presents a review of eigenproblem solution techniques. The 
methods are evaluated and a selection is made of those methods deemed 
necessary for incorporation in the general FEM system that employs modal 
synthesis for substructure reduction. The second section discusses the 
methods used to solve the differential equations of motion. Emphasis is 
placed on the effects that nonlinear structural response and multilevel 
substructured modeling have on the solution processes. The final sec- 
tion briefly discusses mass matrix formulation, damping, and the use ot 
experimental data. 
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4.2 Eicenproblem Sol ut Ion 
4 .2 . 1 Effects of Multi level Snhstructur I ng 

The Introduction of multilevel substructuring into the FEM has 
added significantly to the requirements for a versatile elgenproblem 
solution package in the general FEM system. The importance of accuracy 
In computing substructure modes and the frequency with which the eigen- 
problem must be solved become major considerations in the selection of a 
solution process. Accurate synthesis of the highest level structure 
stiffness and mass is dependent upon the quality of the retained modes 
from each substructure. In a model with multilevel substructuring, the 
stiffness and mass matrices at the highest level are synthesized using 
vibration modes from the lower levels. As each higher level substruc- 
ture Is assembled from the lower levels, repeated vibration analysis ot 
synthesized substructures Is performed. This process will tend to com- 
pound any errors that may exist in the vibration modes. Eventual 
degradation of the highest level stiffness and mass can be avoided by 
maintaining good accuracy in the computed mode shapes at each level. In 
the modal synthesis process, a free vibration analysis is performed on 
each substructure. In nonlinear problems, the vibration analysis ot 
nonlinear substructures is repeated with each update of system proper- 
ties. For large substructures these requirements for frequent vibration 
analyses could become prohibitively expensive unless very efficient 
solution methods are used. 
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In the evaluation of elgenproblem solution methods, one must con- 
sider more than the computational efficiency with which a solution Is 
achieved. Also of prime concern Is the stability of the method for a 
particular problem, given the form of the matrices involved. If the 
solution method Is Inappropriate for the given problem, convergence to 
the wrong solution or divergence can occur. This stability is of par- 
ticular Interest for finite element models using multilevel substruc- 
turing. As the synthesis progresses upward through the various levels 
of substructuring, the nature of the elgenproblem changes as outlined 
below. Thus the analyst must have the capability to select the solution 
method that most closely matches the characteristics of the new problem. 

The changes In the nature of the elgenproblem Include variations in 
the characteristics of the substructure stiffness and mass matrices and 
in the number of eigenpairs needed. For lowest level substructures, the 
mass matrix may be diagonal or banded, positive definite or semt- 
definite, well formed or ill-conditioned. Although the stiffness matrix 
is normally banded. It may be singular, such as when rigid-body modes 
are Included. Further complications arise when static condensation is 
used to eliminate massless degrees of freedom. This may cause a loss of 
bandform for both the stiffness and the mass. At the lowest substruc- 
ture level, usually only the lowest frequencies and mode shapes are 
necessary to progress with the synthesis process. However, depending 
upon the elgenproblem solution method used, this requirement may Involve 
finding all eigenpairs of a reduced system. 
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When a structure's stiffness and mass are assembled from lower 
level substructures, the resulting matrices can take on a special form. 
If Guyan reduction Is used at the lowest level, the structure matrices 
are generally full or widely banded. If modal synthesis Is used, the 
stiffness matrix will be In block diagonal form and the mass matrix will 
be block diagonal and border banded; see Eq. (3.25) and (3.26). It Is 
most desirable to have elgenproblem solution methods that can deal ef- 
fectively with these special forms when the total number of DOF Is 
large. 

No single eigenvalue solution method has been developed which is 
most efficient for all of the expected forms taken by the elgenproblem. 
Within a particular class of problems (e.g. all elgenpalrs are required 
and the system matrices are narrowly banded) superior solution methods 
can be Identified. An adequate selection of these methods must be 
available In a general FEM system, thus permitting the analyst to tailor 
the solution to the particular characteristics of his problem. 

4 .2 .2 Solution Msthflda 

For discussion purposes, eigenproblem solution methods are 
separated Into three categories. These are simultaneous iteration, 
transformation, and polynomial /vector Iteration methods. Solution tech- 
niques within each of these categories rely on one or more of the fun- 
damental properties of elgenpalrs. Some of these properties are ex- 
pressed in equation form as 

1. DG(<j>} s (4.1) 

This form suggests an Iterative solution by assuming a mode 
shape and successively Improving the approximation until con- 
vergence to the true mode shape Is attained. 
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2. M T [K][>] = [A] and (4.2) 

M T [M]M « [I]. (4.3) 

This form suggests that, since the vectors in [<f>] are unique 
within a scalar multiple of themselves, successive transforma- 
tions of DO and CM] To bring them to diagonal form will yield 
the matrices [<|>] and [x]. 

3. det(K - AM) = 0. (4.4) 


This form shows that the roots of the characteristic polyno- 
mial define the eigenvalues of the problem [K]{cj>} = A[M]{cj)}. 


4. 


A (m> < A (m) < A (m+1 ) , < A <m+1) A (m) (4.5) 
1 '• 1 ~ 2 ~ 2 ~ ~ n-m-1 “ n-m 

where A^ and x! m+1 ^ define the eigenvalues of the "m" th and 
the "m+i" th constraint problem for mode "i". This eigenvalue 
separation property states that the characteristic polynomials 
p(A (m) ); m-1,2,...,(n-1) form a Sturm sequence. 


Details regarding these properties may be found In the references [4.2, 
4.25, & 4.27]. 


4. 2 .2.1 Simultaneous Iteration 


Simultaneous Iteration methods Involve successive improvement of a 
set of orthogonal vectors chosen to approximate the mode shapes. Each 
of the methods uses one of the basic forms of inverse Iteration: 


Cx.. f ] - DODO, or (4.6) 

J+1 J 

[K][x j+1 ] = [M][x.]. (4.7) 


In these forms [K] Is the (nxn) structure stiffness, 
structure mass matrix, [x .] Is the (nxq) matrix of 
for iteration "j” and in symbolic form: 

[A] = [L] ^Cm]CL] T , where 

DO = [l]DJ t . 


CM] Is the (nxn) 
assumed mode shapes 


(4.8) 

(4.9) 
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The differences In the methods lie In the processes used to or- 
thogonal ize CXj +1 H and to accelerate convergence. 

Most of the simultaneous Iteration methods reviewed use the stan- 
dard form of the eigenvalue problem, Eq. (4.6). Jennings [4.121 
presented the first such practical method In which the coupling among 
the vectors of [x .1 is evaluated and reduced by an interaction analysis. 
The process involves examining the off-diagonal terms of the product 


[B] = [x.] T [x j+1 ], 


(4.10) 


then decoupling the vectors in [x. + jl and reorthogonal izing prior to the 
next iteration. 

Rutishauser [4.231 introduced a procedure in which the matrix of 
modal vectors is improved by solving the reduced (qxq) elgenproblem 
using the matrix 


ra = Cx J+1 3 T K J+t :. 


(4.11) 


No explicit orthogonal I zat I on is then required. The remaining varia- 
tions of these methods, primarily attributable to Jennings and his col- 
leagues [4.5, 4.7, & 4.131, represent modifications to enhance conver- 
gence or to improve the orthogonal izat ion process. 

Dong [4.8l and Bathe [4 .21 presented independent developments of 
the method most commonly known as subspace iteration. This method uses 
the form of Eq. (4.7) and relies on the solution of a reduced eigen- 
problem similar to that of Rutishauser. This method has a significant 
advantage in that it does not require conversion of the elgenproblem to 
the standard form of Eq. (4,6). This feature affords the potential to 
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solve a wider range of problems (QG and/or DO ill-conditioned) than 
was previously feasible. Several methods have been proposed to ac- 
celerate the convergence of the subspace iteration process [4.1, & 

4.283. These Include shifting, overrelaxation, and the use of Chebyshev 
polynomials. An improved method for selecting starting vectors is also 
available. Although these techniques were Introduced for application to 
subspace iteration, they are applicable to other simultaneous Iteration 
methods. 

4 .2 .2 .2 Transformation Methods 

Transformation methods are based on the orthogonality property of 
eigenvectors. The approach Is to generate a sequence of orthogonal 
transformation matrices, [P.3, of order (nxn), that drive the symmetric 
matrix, [A3, also of order (nxn), to diagonal form. For a case re- 
quiring n k n iterations to achieve convergence, the operation takes the 
form 

CP k D T CP k . 1 ] T ...CP 2 ] T CP,] T [AXP 1 :CP 2 3...CP k ] = M. (4.12) 

The diagonal matrix [A. 3 contains the complete set of eigenvalues. The 
corresponding eigenvectors are computed by 

[cf>3 = [P 1 3[p 2 3 ... [P k 3 . (4.13) 

The eigenpairs are not ordered as In the case of simultaneous iteration. 
Further, since the complete elgensystem Is solved, the resulting 
diagonal matrix and the product in Eq. (4.13) are unique. This unique- 
ness exists only when the transformation matrices are of order (nxn). 
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The generalized Jacobi method Is recognized to be the most ef- 
ficient, truely Iterative, transformation method for solving the general 
el gen problem 

CK]{cf>} - ACM] {<(>}. (4.14) 

Each transformation matrix, CP.], performs a rotation which zeros cor- 
responding off-diagonal terms in DO and CM]. One sweep is completed 
when each off-diagonal term in Q<] and CM] has been individually zeroed. 
Since zeroing one off-diagonal term tends to make a previously zeroed 
term nonzero, multiple sweeps are required to achieve convergence. 

If the eigenproblem can be transformed to standard form: 

CA]{cj>} = A{cj>}, (4.15) 

where CA] is defined by Eq. (4.8), transformation of CA] to tri diagonal 
form can be performed without iteration. Givens Introduced a method of 
plane rotations C4.25] In which, for step "j"» zeros are Introduced in 
row "j” and column "J" without destroying the zeros from previous steps. 
The eigenpairs of the resulting tridiagonal matrix can then be found 
with relative ease. Householder developed an improved approach to 
tridiagonal ization using reflection matrices to perform the transforma- 
tion. Like that of Givens, Householder's method requires (n-2) steps to 
complete the reduction but each step of Householder's method involves 
roughly half as many multiplications. 
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A third method for tridiagonal Ization is available and can be ap- 
plied to the general eigenvalue problem of Eq. (4.14). The generalized 
Lanczos method [4.24H is applicable to problems when QG and [M] are 
symmetric^ equally banded, and at least one Is positive definite. In 
the method, a single orthogonal transformation matrix is built one 
column at a time. With the evaluation of each new column, {r ; }, reor- 
thogona I ization to all previous columns, {r^ to (r ; _^}. Is required to 
ensure stability. 

Solution for the eigenpairs of a tridiagonal matrix is effectively 
attained by using QR iteration to find the eigenvalues and Inverse vec- 
tor iteration to find the eigenvectors. A discussion of Inverse vector 
iteration follows in the next section. 

QR Iteration C4.21 is simply a factorization of the tridiagonal 
matrix, CS], Into the product of an orthogonal matrix, CQD» and an upper 
triangular matrix, CR], 

CS] !S [Q][R] (4.16) 
Pre- and post-multiplying Eq. (4.16) by [Q3 and EQ] respectively gives 

naTsiM = crim (4.1-/) 

where the product on the right-hand side is the diagonal matrix of 
eigenvalues. Jacobi rotation matrices are often used to reduce [SU to 
uppertriangul ar form of DO, with []Q3 being the product of those rota- 
tion matrices. A shifting strategy, discussed below, can be used with 
QR iteration to accelerate convergence. 
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4.2 .2.3 Polynomial /Vector Iteration 


Polynomial Iteration techniques involve an iterative solution of 
the characteristic polynomial 


p(A) = det([K] - A[M]) = 0. 


(4.18) 


In practice the polynomial Is never explicitly evaluated. Instead an 
estimate of the exact eigenvalue Is made and then the determinant is 
evaluated using Gauss or Choleski factorization. From this evaluation, 
the estimate of the eigenvalue is improved and the process Is repeated. 

Two popular procedures aid in this evaluation; accelerated secant 
iteration and bisection. Accelerated secant Iteration Is a super I inear 
Interpolation method from which an improved estimate of the eigenvalue, 
A k+1 » Is determined using two previous estimates, A^ and A^. The im- 
proved estimate is 


A 


k+1 


a 


fpa k )(X k ~ X k-1 

f ( y - p(A k-i ) 


> 


i 


(4.19) 


where a is an acceleration constant. The polynomial is evaluated at 


k+1 


and If required, another iteration is performed. Bisection is an 


iterative process In which the interval between the previously deter- 
mined upper and lower bounds of A^ Is halved and the sign of the polyno- 
mial at that midpoint evaluated. On the basis of the evaluation, the 
new interval containing A Is Identified and the process repeated until 

K 

convergence is achieved. 
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Vector Iteration Is a special version of simultaneous iteration in 
which only one eigenpair Is computed at a time. The iteration sequence 
takes the form: 


DCKx } = CM]{x.}. (4.20) 

J +1 J 

If convergence of the eigenvector is gained In Iteration "k", the eigen- 
value can be computed from Rayleigh's quotient: 


{x k } T M(x k } 

{x k } T [M]{x k } 


(4.21) 


The above form of vector iteration, termed Inverse Iteration, will Ini- 
tially converge to the most dominant eigenpair (smallest eigenvalue). 

To force convergence to a subdominant eigenpair, the Influence of 
the previously computed eigenvectors must be eliminated. If an estimate 
of the eigenvalue Is known, shifting can be used to make the unknown 
eigenpair dominant. Letting y be the value of the shift, a modified 
el gen problem Is defined 


D< - yMH(cj)} = vOGW 


(4.22) 


where v = X- y. inverse Iteration applied to this problem produces the 
eigenpair closest to the shift y. When a good estimate to the desired 
eigenvalue is unavailable, shifting Is not effective. The Iteration 
vector must be chosen to be orthogonal to the previously computed eigen- 
vectors. When eigenvalues are multiple or clustered. It Is necessary to 
reorthogorial ize the iteration vector to avoid convergence to mode shapes 
which have already been computed. 
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Effective eigenproblem solution techniques using a combination of 
polynomial iteration and vector Iteration have been developed for 
limited applications, in Gupta’s method [4.103 the individual intervals 
containing each eigenvalue are defined using Sturm sequence properties 
and bisection. A shift is made to the middle of each interval and then 
the eigenpalr is determined by inverse iteration and Rayleigh's 
quotient. The upper bound on the interval for A. is used as the lower 
bound for A. + ^. Bathe and Wilson presented a method called determinant 
search [4.2J In which accelerated secant iteration Is used to accurately 
determine the eigenvalue, inverse iteration with shifts is then used to 
find the eigenvector. If the eigenvalue is accurately evaluated, a 

shift to that value will result in the eigenvector being found In no 
more than two vector Iterations. An advantage of these combined tech- 
niques Is that when eigenvalues are not multiple or clustered, each 
eigenpalr Is Independently evaluated with no need for orthogonal izat I on. 

Robinson and Harris [4.223 Introduced a Newton-Raphson Iterative 
approach to vector iteration. For Iteration ’’i”, the process involves 
choosing a vector Increment, {Ax.}, and a frequency increment, AA., 
which will eliminate the residual vector 

frj} » [K3{X;} - A.[M3 {x.}, (4.23) 

and thus force convergence to a true eigenpalr. As an additional side 
condition, the vector increment, {Ax.}, must be orthogonal to the ap- 
proximate eigenvector, {x }, with respect to [M3, i.e. 

{X.} T [M3{AX.} = 0. 


(4.24) 



For DO and DO of order (nxn), a set of (n+1) simultaneous equations is 
then solved for {Ax.} and AA.o The process Is repeated until some con- 
vergence criterion Is satisfied. Upon convergence, estimates of a dif- 
ferent elgenpalr are taken and the process is restarted. Initial 
estimates to the eigenpairs are taken from another source, such as sub- 
space iteration. 

The above method requires a complete trt angulation of an (n+1) 
order set of linear equations for each iteration. Lee and Robinson 
[4.17] presented an alternative method which requtres only one 
triangulation for each elgenpalr. Although the convergence rate Is less 
than for the Robi nson-Harris method, the overall computational efficien- 
cy Is improved. When eigenvalues are multiple or clustered, all eigeri- 
palrs in the group are found simultaneously by a method using Lagrange 
multipliers and the stationary property of the Rayleigh quotient. Il- 
lustrative examples show the Lee-Robinson method to require between 36? 
and 60? as many operations as needed by subspace Iteration. Comparison 
of operation counts Indicates that the method Is also more efficient 
than determinant search. 

4.2 .3 Eval uation of Eiqenprobl em Sol ution Methods 

As with the evaluation of modal synthesis, it is useful to 
establish criteria for the evaluation of eigenproblem solution methods. 
Although the categories listed below are similar to those previously 
discussed, they differ in the application to eigenproblems. 
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1 . Computational £ff,or.t 

The computational effort required by a solution algorithm 
is most effectively measured by the number of arithmetic 
operations performed. For elgenprobiem solutions, this number 
is dependent upon the order and bandwidth of the matrices and 
the number of eigenpairs computed. For large problems. In- 
core solutions are not always possible. Therefore, efficiency 
is further measured by the capability to solve the eigen- 
problem with the use of out-of-core data storage. Controlling 
parameters are the number of data blocks required to simul- 
taneously reside In core and the frequency with which these 
blocks must be replaced. 

2. Appl I cab i I ity ±e Genera I Prob I ems, 

Generality of eigenproblem solution methods Implies the 
ability to solve problems of the form given by Eq. (4.14). 
Often the transformation to standard form, Eq. (4.15), Is 
ineffective because of an ill-conditioned mass matrix. Such a 
situation Is common with a lumped mass formulation in which 
rotational DOF are assigned zero inertia. This limitation is 
avoided when an eigenproblem solution method Is available that 
can solve the general problem. Generality further applies to 
dealing with the various forms taken by the stiffness and mass 
matrices. As shown in Chapter 3, the dynamic reduction 
process can lead to coefficient matrices that are diagonal, 
banded, full, b I ock-di agonal, or block- diagonal and border 
banded. 

3. at Reanai y,sis 

The need for reanalysis of an eigenproblem arises in two 
Instances. First, the analysis may need to be continued with 
a tighter convergence limit to attain greater precision in the 
frequencies and mode shapes. Secondly, as the physical 
characteristics of the model change during nonlinear response, 
the vibration characteristics also change. V/hen a previous 
eigensol ution Is available, significant cost savings are 
realized if the method takes advantage of this information in 
reana lysis. 

4. A-C-C.ura.cy and Stab 1 1 1 tv 

As with modal synthesis, stability of the operations and 
accuracy of results of the eigenproblem solution methods are 
vital. Stability problems tend to be caused by mode shape 
vectors that are not orthogonal while accuracy is controlled 
by the convergence criteria. 



The superior eigenproblem solution methods from each of the three 
categories are Identified below. The use of these methods and their In- 
terdependencies will be discussed In the next section. 

Subspace Iteration emerges as the superior process In the category 
of simultaneous Interation. The primary advantage over the Jennings and 
Rutshauser approaches Is Its applicability to the general problem form. 
The potentially ill-conditioned transformation to standard form Is 
avoided, thus making the method generally more stable. Other advan- 
tages, which are common to all methods in this category. Include the 
ability to use previous modal Information In reanalysis and the ease 
with which acceleration schemes are incorporated. 

The most efficient solution process of the transformation methods 
category is that employing the Householder transformation to tr I diagonal 
form followed by QR Iteration (the HQRI method). This method is effec- 
tive only when the eigenproblem can be written in standard form and when 
alll eigenpairs must be calculated. Although the method does not take 
advantage of the bandform of the equations, this is not necessarily a 
drawback since the conversion to standard form often results in a fully 
populated coefficient matrix. Such is the case when a consistent or a 
synthesized mass matrix is used. 

When the above conversion cannot be made, general I zed Jacobi Itera- 
tion makes an attractive alternative to the HQRI process. Although 
Jacobi iteration is not as computationally efficient as other transfor- 
mation methods, it has two distinct advantages. First, as the off- 
diagonal terms become smaller, the process becomes more efficient, i.e. 
convergence Is more rapid on systems that are almost diagonal. 



Secondly, Jacobi Iteration Is the only transformation method which ef- 
fectively utilizes a previous solution in reanalysis. Rather than 
starting the process with a standard rotation matrix to zero an off- 
diagonal element, the previous mode shape set Is used to rotate the 
coefficient matrices close to diagonal form. 

In general, as the order of the elgenproblem and/or its bandwidth 
increases, polynomial /vector Iteration methods become less economical 
because of their reliance on computing the determinant in Eq. (4.4). An 
exception is the method presented by Lee and Robinson [4.17]. In this 
method only one complete triangulation is performed for each efgenpair 
that is evaluated. The required starting solution, in the form of a 
prior modal analysis, can be used to advantage in reana lysis. 

As mentioned previously, operation counts provide the best evidence 
regarding the efficiency of a solution algorithm. For comparison pur- 
poses, Table 4.1 lists the operation counts for the methods discussed 
above. It is assumed that no computational penalty is paid by Im- 
plementing the methods with hypermatrix data structures. 



Table 4.1 


OPERATION COUNTS FOR E IGENPROBLEM SOLUTION 


Method 

Number of operations required F4.2, 4.17] 

Subspace Iteration!^ 

(2) 

HQRI. . . 

(2 3) 

Generalized Jacobi Iteration. ", 

Lee-Robir.son Vector Iteration 

. . nm^ + nm(4 + 2p) + 4np + 20nq(2m + q + 1.5) 

. . 0.67n^ + 10. 5n^ + pn^ + 9pn 

, 3 2 

. . 3n + 6n 

2 

. . 0.5np(m + 5m + 2) + nT(7m + 6) 

Notes 

Notation 

(1) Assumes 10 iterations required 

n = order of fK] and FM] 

(2) Assumes fully populated FK] and FM] 

m = half-bandwidth of FK] and FM] 

(3) Total for one sweep only 

p = number of eigenpairs computed 


q = mi n(2p, p + 8) 


T = number of iterations required 
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4.2.4 Choice of Elgenproblem Solution Methods 

It has been demonstrated that a wide variety of el genprob terns can 
arise In dynamic analysis by the FEM. As mentioned earlier, no single 
method can be effectively applied to the solution of each of these 
problems. Therefore, a selection of methods must be available that per- 
mits the analyst to choose the one most suited to his problem. The 
three controlling variables in making the choice are the number of DOF 
in the model, the bandwidth of the matrices, and the number of eigen- 
pairs to be computed. 

In general solution methods from the three catagorfes are most ef- 
ficient when applied to a particular class of problem. When all eigen- 
pairs are required of large, fully populated matrices, a transformation 
method should be used. When only a few frequencies and mode shapes are 
needed and the equations are narrowly banded, polynomial /vector Itera- 
tion methods are most efficient. Simuntaneous iteration exhibits Its 
superiority in the range of problems between the two mentioned above. 

To fulfill the selection requirements, it is proposed that the 
superior methods identified in the preceding section be Incorporated in- 
to the general purpose FEM system. These methods ares subspace itera- 
tion, HQRI, generalized Jacobi Iteration, and the Lee-Robinson vector 
iteration. Conveniently, several of these methods can be used In con- 
cert. That Is, a starting solution for Lee-Robinson vector iteration 
may be obtained from subspace iteration. Also, subspace iteration ef- 
fectively utilizes generalized Jacobi iteration in the solution of the 
reduced elgenproblem. When the transformation to standard form, re- 
quired by the HQRI approach, is not effective, generalized Jacobi itera- 



tion becomes the next best alternative. 


Although operation counts provide useful Information in choosing a 
solution method, other factors must be included. The availability of a 
starting or Initial solution may significantly increase the efficiency 
of some methods. Also, the number of iterations required by some 
methods Is highly dependent on the convergence tolerance that is 
established. These factors force selection of the appropriate eigen- 
problem solution method to be made more on the basis of the analyst’s 
experience and judgement rather than on published operation counts. 



4.3 Sol utlon of the Equatl ons of Motion 


4.3.1 Introduction. 

Consider the form of the differential equation of motion for a 
damped finite element system; 

[M(t)]{u} + [C(t)]{u} + [K(t)]{u} = (P(t)}. (4.25) 

For many problems the coefficient matrices of Eq. (4.25) are constant 
with respect to time, resulting In a set of coupled linear differential 
equations. When material and/or geometric non II near I ties do arise, they 
are most often I Imlted to the stiffness matrix, CK(t)3. The mass of the 
structure seldom varies with time and the damping characteristics 
generally are not sufficiently well understood to warrant any time 
dependent change In the modal damping ratios. 

When multilevel substructured models are used and the coefficient 
matrices are derived by modal synthesis, the mass and damping matrices 
can also become time dependent. These dependencies enter the syn- 
thesized mass through the coordinate transformations performed on each 
condensed nonlinear substructure. Likewise, the relationship between 
damping and the structure mass (see Section 4.4) causes a corresponding 
variation In the damping matrix. Therefore, In Its most general form 
Eq. (4.25) has all three of Its coefficient matrices varying with time. 

The transient response of a substructured model is described by the 
displacements {u}, velocities {u}, and accelerations (u } for the DOF In 
the highest level structure. Two general strategies are available for 
obtaining this solution. They are mode superposition and time-history 
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Integration. Both strategies are fundamental to the dynamic analysis of 
structures and as such are familiar to most analysts. 

The selection of a particular solution strategy is highly problem 
dependent. The decision variables include the degree of nonlinearity 
(If any), the number of modes participating in the response, the length 
of the time interval over which the response must be evaluated, and the 
nature of the transient load. A discussion of these variables and of 
the effects of multilevel substructured modeling on the solution 
procedures is presented In the following sections. 

4.3.2 Mode Su pe r posl t lon 

Mode superposition is an economical approach to solving the equa- 
tions of motion when the structural response is linear and Is limited to 
the lower modes of the frequency spectrum. For linear systems, the 
basic approach Is to uncouple the equations by a transformation to modal 
coordinates. The resulting set of equations describe the dynamic 
response of the structure in each mode when excited by the corresponding 
modal component of the applied load. These modal equations are solved 
Individually by a direct Integration method and then the total response 
is computed as the sum of the modal responses. 

A computational advantage with modal superposition over time- 
history Integration methods (see Section 4.3.3) is realized when the 
dominant portion of the structural response is contained in the lower 
modes. This characteristic Is evident when the loading function has 
very little high frequency content. Examples of such loading types 
might include earthquake excitation and mechanical vibrations. The 
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number of modes that participate In the response and thus must be In- 
cluded In the modal analysis Is easily estimated by examining the 
loading function after It has been transformed to modal coordinates. 
Loadings that typically have more high frequency content Include blast 
or shock loading and Impact. These later examples are more economically 
.analysed by time-history Integration of Eq. (4.25). 

The extension of mode superposition to nonlinear systems has 
received only limited attention In the literature. Morris C4.18H has 
applied the concept to cable systems and framed structures. His results 
are not particularly encouraging. The method performed poorly when ap- 
plied to structures experiencing more than Just mild non 1 1 near I ties. 
Nickel I C4.213 presented a more complete development of nonlinear mode 
superposition. Although mode superposition performed very well In his 
example problems* It showed no particular numerical advantage over 
direct Integration schemes. Nickel I does note that some physical In- 
sight Into the behavior of the structure Is gained by observing the 
changes in the modal spectrum as deformation proceeds. 

The procedure for solving Eq. (4.25) by mode superposition remains 
unchanged when multilevel substructured models are Introduced. The 
transformations performed on the coefficient matrices are Independent of 
the methods used to build those matrices. Although both geometric and 
modal DOF from the condensed substructures may exist In the structure 
stiffness* damping* and mass matrices* no complication is expected in 
making further transformations to uncouple the equations of motion. In 
fact* the transformations performed in mode superposition are quite 
similar to those used In modal synthesis. However* no studies have been 
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located that document the effects of substructuring on the mode super- 
position solution strategy. 

4.3.3 T I trie- Hi story I ntearation 

When the dynamic loading on a structure excites many modes ot 
vibration, time-history Integration, also known as temporal or direct 
Integration, Is often a more economical solution procedure than mode 
superposition. The Integration operators work directly on the coupled 
equations of motion without regard for the modal content of the 
response. 

Operators for time-history Integration vary In the relationships 
between the known and the unknown solution variables. These variables 
are the displacements, velocities, and accelerations of the structure 
DOF at various points In time throughout the response period. The as- 
sumed relationship between the known and unknown solution variables 
defines whether the operator Is explictt or Implicit. In explicit 
methods [4.2 J, the displacements at the end of the tnterval, {u ++A+ }, 
are based on the equilibrium conditions at the start of the interval: 
Eq. (4.25) evaluated at time t. After a starting solution has been 
established. It Is possible to progress through the solution without 
solving a set of simultaneous equations. This feature of explicit 
methods presents an advantage In computational efficiency over the im- 
plicit methods discussed below. However, since explicit methods are ex- 
trapolatory In terms of the satisfying equilibrium, they are only con- 
ditionally stable [4.15] at best. A stable integration operator Is one 
for which the error Introduced In each time step Is bounded for the 
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chosen step size and thus the total error does not grow without bound as 
the solution progresses. A conditionally stable operator is one which 
remains stable for only a limited time step size. The central dif- 
ference method [4. 9], when used with a diagonal mass matrix, is con- 
sidered to be an accurate and economical explicit solution method 
■ [4.163. 

Implicit operators rely on the solution of Eq. (4.25) at time 
t + At to obtain the solution {u }. Two significant consequences 
result. First, a set of linear equations must be solved at each incre- 
ment In time considered In the solution. Secondly, ft is possible to 
develop Implicit operators that are unconditionally stable (stable for 
any size of time step At). The two most popular implicit operators are 
the Newmark-8 method [4.193 and the Wilson-9 method [4.263. 

The Implicit operators which exhibit unconditional stability when 
used on linear problems may not always maintain this feature when 
general nonlinear response occurs. Proofs of unconditional stability of 
the Newmark-8 operator, when applied to nonlinear problems, have been 
presented [4.4 & 4.73 but contradictory evidence has also been 
documented [4.203. The question of unconditional stability Is not 
necessarily a drawback to the use of Implicit (or even explicit) 
operators. For nonlinear analysis relatively small time steps (which 
lead to small displacement changes) must be taken to accurately detect 
and incorporate nonlinear behavior. As a result, frequent equilibrium 
iterations with updates of the coefficient matrices will be required to 
maintain accuracy, regardless of which Integration operator Is employed. 
This process can succeed in arresting any instability that may develop 
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from using a particular integration operator [4.203. 

Explicit integration of the equations of motion has proven to be 
effective for nonlinear problems that do not involve substructured 
models. Since the frequent equilibrium iterations prevent instability 
in the solution, an implicit operator that may be unconditionally stable 
will not show any advantage over explicit operators. Further, It is not 
always necessary to use the structure stiffness matrix with explicit 
operators to compute the internal nodal forces. Computation of the 
elastic force vector can be performed directly by integration of the 
following expression for each element in the structure [4.33: 

(F e > = X[L3 T | v [B3 T {o}dV, (4.26) 

in which {F q } = structure elastic force vector, 

[L3 = connectivity matrix for the element, 

[B3 = strain-displacement relation, 

{a} = element stress vector, 

V = volume of the element, and 
n = number of elements in the structure. 

In substructured models the element stress vector is not available 
for condensed lower level substructures. As a result, computation of 
the elastic force vector must incorporate at least a portion of the 
structure stiffness matrix. The resulting elastic forces are 

tF e> = {F e>H + K>c {u> c <4 - 27> 

where (F } = structure elastic force vector, 
e ' 

(F }., = internal forces for uncondensed elements at the 

e H 
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highest level, computed by Eq. (4.26), 

[K] = stiffness matrix for the adjoining condensed 

c 

substructures, and 

{u} = condensed DOF displacement vector. 

c 

Computation of elastic forces by Eq. (4.27) requires that the stiffness 
matrix for the condensed substructures be brought Into memory during 
each time step. It Is yet unknown whether this requirement will reduce 
the computational effectiveness of explicit Integration of the nonlinear 
equations of motion. 

Since I mp I id t methods require equation solving with the full 
structure stiffness matrix, no comp! I cat ions in the algorithms due to 
the use of substructured models are anticipated. 

4.4 Additional Computational Cons! derations 

Several minor computational details regarding the dynamic analysis 
of structures are complicated by the use of substructured finite element 
models. Details considered in this study include mass matrix formula™ 
tlon, damptng, and the use of experimental data. 

Two approaches are used to assemble the mass matrix of a substruc- 
ture containing only simple elements C4.63« The most elementary ap- 
proach is the lumped mass formulation. With this approach, the mass of 
each element is lumped at the nodal points. Upon assembly, a diagonal 
substructure mass matrix is achieved. This form has many computational 
advantages as previously noted. The alternate approach is the consis- 
tent mass formulation. In this approach, accelerations within the ele- 
ment are Interpolated by the same shape functions used to describe the 
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element displacements. Thus, the assembled substructure mass matrix has 
the same banded nature as the substructure stiffness. Regardless of 
initial formulation, when multilevel substructuring is introduced, the 
resulting structure mass matrix Is of banded form. It is anticipated 
that as the number of levels of substructuring increases, the dif- 
ferences between the lumped and consistent formulations will diminish. 
The reduction process used, through the various transformations per- 
formed, Is expected to mask the original mass matrix formulation and 
thus make the two formulations indistinguishable. 

Structural damping is a phenomenon which Is not easily modelled by 
analysts. Usually, damping can only be evaluated by experimental 
testing thus making Its specification for a general finite element model 
difficult. However, there are procedures designed to include the ef- 
fects of damping In a computationally efficient manner £4.6]. Normally, 
modal damping ratios are known (or selected) and used directly In com- 
puting structural response In modal coordinates. If It Is desirable to 
have an explicit damping matrix in geometric coordinates, one can be 
formulated which Is proportional to both the stiffness and mass 
matrices. In this form, termed Rayleigh damping, the system vibration 
mode shapes are orthogonal to the damping as well as the mass and stiff- 
ness matrices. In substructured models, application of damping is 
delayed until the equations of motion for the highest level structure 
are formulated. The methods of Guyan reduction and modal synthesis do 
not provide for damping at the substructure level. 
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Pre-existing experimental data regarding the free vibration charac- 
teristics of a substructure could be useful for reducing the com- 
putational expenses of a finite element analysis. It is a simple matter 
to adapt input translators to accept this data and thus eliminate un- 
necessary vibration analyses. An important consideration in using ex- 
perimental data is the consistency in boundary restraint between the 
test fixture and the finite element model. 
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CHAPTER 5 


MATRIX FORM OF NONLINEAR CONTINUUM MECHANICS 


5.1 I ntroduct ion 

The dynamic solution procedures described in the preceeding chap- 
ters are independent of the continuum theory used to derive the govern- 
ing nonlinear equations of motion. This chapter focuses on the essen- 
tial aspects of nonlinear continuum theories. This information is in- 
cluded in the report for three reasons: (1) the formulations are not 
readily available in matrix notation familiar to finite element workers, 
(2) the exact transformations in matrix form will be required during the 
implementation phase, and (3) a relative assessment of efficiency of the 
formulations is desired. The finite element formulations are developed 
and contrasted as to their suitability for substructured dynamic analy- 
sis in the next chapter. 

Initial work on a matrix formulation for nonlinear continuum theory 
was conducted by Nayak [5.1]. Zienkiewicz [5.2] has been the primary 
advocate of this approach which has been adopted by most researchers at 
Swansea University. This chapter summarizes their work and presents an 
extension of the approach to include an Updated Lagrangian formulation. 

Nonlinear behavior includes that due to constitutive relationships 
(material nonlinearity) and strain-displacement relationships (geometric 
nonlinearity). To achieve maximum generality in the discussion, struc- 
tures are assumed to exhibit both types of nonlinearity simultaneous- 
ly. The discussion of material nonlinearity is directed towards metal 
plasticity at both infinitesimal and finite strain magnitudes. While it 
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is recognized that strain-rate dependent behavior may be of importance 
in dynamic analysis, these effects are not considered in this discus- 
sion. Few procedures in the literature include strain-rate dependen- 
cies, creep and thermal dependencies of material properties. Most 
often, rate effects are incorporated into the basic concepts of incre- 
mental plasticity and thus do not affect the comparison of nonlinear 
formu I at ions. 

The dynamic analysis of structures in which nonlinearity is limited 
to the material stress-stra i n behavior presents no complications in for- 
mulating the equations of motion. This follows as all quantities are 
referred to the initial configuration of the structure, with displace- 
ments and rotations assumed to be infinitesimal. The challenge is de- 
veloping realistic constitutive models for the material behavior under 
cyclic loading. Geometric nonlinear effects are important when the 
structure undergoes large rotations and/or finite strain magnitudes. 
This occurs, for example, in the following situations: crashworthiness 
of aircraft structures and components, metal forming processes, necking 
of structural components subjected to tensile overloads, geometry 
changes under service loading that affect performance (e.g., turbine 
blades) and the localized deformation of material in the vicinity of 
stress concentrations — especially notches and cracks. 

The governing equations of nonlinear solid mechanics are summarized 
Jn the following sections in terms of Cartesian reference axes. The 
notation adopted here is a modification of that first used by Bathe 
[5.3!. Each quantity retains the same symbol throughout the deformation 
process. Superscripts and subscripts are used to indicate the config- 
uration of the structure in which the variable occurs and to which con- 



figuration it is referenced, respectively. A single left subscript 
implies an increment of the variable over time At with the subscript 
indicating the reference configuration. A summary of the notation is 
presented in Section 5.9. 

A! I equations are derived using matrix, rather than tensor, nota- 
tion. The theory of finite deformation is elegantly expressed in tensor 
notation. However, all tensor operations require transformation to 
equivalent matrix form for computer implementation. This, combined with 
the exclusive use of matrix notation in finite element research, encour- 
ages the adoption of matrices at the start. The reader already familiar 
with tensor notation will readily recognize the appearance of symmetric 
strain and stress "matrices." 

The following sections contain discussions of coordinate systems, 
strain-displacement relations, stress definitions, and virtual displace- 
ments for finite deformations. Particular emphasis is placed on the 
finite deformation strain and stress rates that are required for virtual 
work arguments. Two nonlinear formulations arise in the discussion. 
These are termed Total Lagrangian (T. L.) and Updated Lagrangian (U. 
L.). In the T. L. formulation all quantities are referenced to the con- 
figuration at time t=0. The difficulty arises with constitutive 
relations for the material which are naturally expressed in terms of the 
deformed configuration, e.g., true stress and true strain. The correct 
transformations of strain, stress, strain rates, stress rates, and 
Incremental constitutive relationships between the initial and 
instantaneous configurations are derived in this chapter. The U. L. 
formulation refers all quantities at time t+At to the configuration at 
time t. Numerous attempts have been made to derive an U. L. formulation 
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for finite element analysis. Unfortunately, very few have followed the 
rigorous equations of continuum mechanics; most formulations are derived 
by assuming certain nonlinear terms are small during the motion from t 
to t+At. This has led to considerable misunderstanding in the liter- 
ature. The exact U. L. formulation is rigorously derived in this 
chapter,. Only then are several commonly used simplifying assumptions 
described. 

An extensive discussion of constitutive models for finite defor- 
mation elasto-plasticity is included in this chapter. The various 
strain and stress rate definitions and their transformations are 
presented in matrix notation. This material will enable the realistic 
modeling of large deformation effects that occur near impact zones. It 
is anticipated that the use of substructur i ng prior to a transient anal- 
ysis will show a significant computational advantage in such problems. 

The chapter concludes with a qualitative assessment of the computa- 
tion efficiency of the two major formulations. 

5 . 2 Coord i nate Systems and Transformations 

The configuration of the body is considered at three times; namely, 
0, t, and t+At. For static analysis, the parameter t may be associated 
with a loading state rather than time. Motion of the body is described 
with respect to fixed Cartesian axes having unit vectors {i}. The 
position vector of any material point at the three times is denoted 
°{x> , {x}, and ^" + ^{x). Coordinates °{x} remain fixed for a particular 

material point and thus are convected in the present usage. No 
subscripts are necessary for the position vectors 
correspond i ng differentials. 


and their 
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Let °{dx} represent 

the 

components of 

an initial 

1 i ne 

segment, °ds, which deforms 

i nto 

line segment ^ds 

at time t. 

Compo- 

nents ^{dx} then are given 

by 




+ {dx} = *[J] °{dx} 




(5.1) 


in which ^"[J] 
o 


is the deformation Jacobian with terms defined by 


t 

o 




(5.2) 


The subscript outside the ( ) explicitly indicates the reference 
configuration for differentiation. The deformation Jacobian fully 
character i zes motion in the differential neighborhood of a point that is 
displaced from °P to "*"p. Geometric considerations [5.1! show that 
differential volume and area changes may be expressed as 


+ dV = + I J | °dV (5.3) 

o 1 1 

f {dA} = + I J ! f [J]" T °{dA} (5.4) 

o 1 ' o 

in which °dV is a differential volume and °{dA} is the vector of 
components for a differential area °dA, both at t=0. Vertical bars are 
used to indicate a determinant. 

The position vector ^”{x} is written in terms of a displacement 
vector f {u} with components referred to base vectors {i} as 

f fx} = °{x) + + {u} (5.5) 


The deformation Jacobian may then be written in the form 



[I 1 



(5.6) 
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where [j] is termed the displacement Jacobian and 
o 


matrix. Terms of [j 
o J 


t . ,t . 

o J I ,J o l,J 


are defined by 


[ I ] is an identi ty 


(5.7) 


Equations (5.1) - (5.7) provide the basis of stra i n-stress 
transformat ions required in the T. L. formulation. 

The finite magnitude changes in volume and area of differential 
elements that occur during the motion from t to t+At may be measured 
with respect to the configuration at time 0 or t. First, the config- 
uration az t+At is expressed in terms of the configuration at t and the 
displacement increment (deformation plus rotation). No restrictions are 
imposed on the magnitude of the displacement increment during At. In 
analogy with (5.1) an expression of the form 


t+At , , t+At t. 

{dx} = + [J] {dx} 


t+At , 


(5.8) 


is sought in which ^ [ J ] is the increment of the deformation Jacobian 
over At but measured with respect to the configuration at time t. The 
terms of the Jacobian increment are given by 


t+At . ,t+At . 

, J . . = . ( x. .) 
t ij + i,j 


Let {Au} represent the incremental displacement vector: 
1 ' 1 " 1 ° ~ I 

| A u, Av, Awl during motion from t to t+At. Then 


(5.9) 


+ {Au} T = 


t+At t t 

+‘V = t' *i,j> + t ( 4 1 ,j 


(5.10) 


More simply, (5.10) may be written 


t+At t+At 

I r j i = m + iivi 

t t 


(5.11) 
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in which *" +A |[v] represents the mean displacement gradient with respect 
to the current (t) configuration (or equivalently the mean velocity 
gradient if divided by At). Increments of differentia! volumes and 
areas from t to t+At may be similarly expressed in analogy with (5.3) 
and (5.4) as 


t + 4t dv . fW (j| 
t+At i 


t 


t+At 


{dA} 


dV 


i t+At -T t. 

J[ + IJ] {dA} 


(5.12) 

(5.13) 


Equations (5.8) - (5.13) provide the basis for transformations required 
in the U. I. formulation as first presented by Yaghma? and Popov 
[5.4]. The U. L. formulation is generally attributed to Bathe [5.3, 
5.5] due to the numerous papers he has published on the formulation. 

Variables for the configuration at t+At may also be related to that 
at t=0 through a sequence of transformations — first from t+At to t, 
then from t to t=0. Consider first that at t+At, 


t+At. t+At o 

{dx} = q [J] {dx} 

Subtraction of (5.1) from the above expression yields 


(5.14) 


t+At 


{dx} - {dx} = [AJ] {dx} = [ T a [J3 - ' [ J : 
o o 


{dx} (5.15) 


in which [AJ] is the incremental deformation Jacobian but referred to 
o 

the configuration at t=0. In ( 5 . 8 ), ^ [ J ] is the incremental defor- 

mation Jacobian referred to the configuration at time t. Using the 
incremental displacement vector {Au} over the time increment At, the 

terms of [AJ] are 
o 


(AJ. . ) = ( + x. + + Au. ) . . - ( + x. .) 

o ij o i i i»J o l,J 


(5.16) 
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or, more simply 

(AJ . .) = + (Au) . . (5.17) 

o ij o l,J 

Using the chain rule, the derivatives in (5.17) with respect to { x } are 

-J* 

transformed to those with respect to {x}. The operations readily show 
that the result may be written as 


t+At t 

[AJ] = . [ v] [J] 

o to 


(5.18) 


Substitution of (5.18) into (5.15) shows that 


t+At 


{dx} = [ [I ] + ++A !lv] ] t [J] °{dx} 


(5.19) 


and since {dx} = [J] {dx}, the above expression simplifies to 


t+At, , t+At t f . 
{dx} = + [J] {dx} 


(5.20) 


which is the same result obtained in (5.8). This procedure illustrates 
the multiplicative decomposition of deformation that parallels simple 
coordinate transformation (which must be the case when the displacement 
gradients represent a rigid rotation). 

Consider now the imposition of virtual, rather than finite, dis- 
placement increments on the configuration at time t. Correspond i ng 
virtual changes in differential areas and volumes are desired. Given 
the differential volumes and areas at time t, the expressions previously 
derived for finite changes during At are re-examined for the limiting 
case as At * 0, i.e., the virtual displacement occurs during a vanish- 
ing increment of the time (or loading) parameter. introducing the first 
variational operator, <$, the variation of line segments is given by 


5 f {dx} = *[6J] °{dx} 


(5.21) 
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in which terms of the variation of the deformation Jacobian appearing 
above are defined by 

f <$ J . . = ( + 6x. .) = ( + 6u . .) (5.22) 

o ij o i,j o i,j 

since the variation of the initial coordinates, °{x}, is zero under the 
imposed virtual displacement. If the virtual displacements, "*"{5u}, are 
interpreted as occurring over a time dt, the derivatives in (5.22) may 
be considered virtual velocity gradients. Using the chain rule, the 
virtual displacement gradients are referred to the configuration at time 
t with the resulting form 


+ I<$J] = T E <S v ] + [J1 (5.23) 

o to 

in which the variational form of the displacement gradient matrix 
defined in (5.17) and (5.18) is used. Substitution of (5.23) into 
(5.21) provides the variation of the deformation Jacobian with respect 
to the current configuration in analogy with (5.8) 


5 + {dx} = *{5v! *(J3 °{dx} 


(5.24) 


but ^[J] °{dx> is simply + {dx}, thus the above simplifies to 


S + {dx} = ^E5vl + {dx> 


(5.25) 


which implies that 
JjSJJ = Jlfivl 


(5.26) 


Using Eufer ? s theorem for homogeneous functions, it is shown that the 
variation of the deformation Jacobian determinant is given by 


ol«j| = T r(J[«vl) r o U\ 


(5.27) 
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in which Tr( ) denotes the trace of the matrix. It is customary to 
express the virtual displacement gradient matrix, ^[6v], as the sum of a 
symmetric deformation matrix, |[6e], and an anti -symmetric spin matrix, 
^[<5w], as 

Ttfivl = tlfiel + T(6w] (5.28) 

t t t 

in which 

f f lSe] j[ [ 5 v I + [6v! T ] (5.29) 

and 

JlSwI * j *[ [ 5 v I - [6 v] T ] (5.30) 


Both the deformation and spin components are linear in the virtual dis- 
placement derivatives. Biot 1 5.6] has demonstrated that ^[6w] 

1 ” 

represents the rotation of rigid-convected axes attached to point p 

during the small (infinitesimal) virtual displacement. Similarly, 

|[<5e] is shown to constitute the pure deformation caused by the virtual 

displacements. A null deformation matrix is therefore a necessary and 

sufficient condition for a rigid virtual displacement. The 6 x 1 vector 

*}- 

representation of the symmetric deformation matrix is denoted ^{6e} 
(shear terms are doubled to form the vector representation) . 

Using the above procedure, the following relationships are derived 
for the U.‘ L. formulation. 



(5.32) 
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These relationships are analogous to (5.21) and (5.23). The virtual 
displacement gradients at time t+At may again be expressed as the sum of 

++A + 

a symmetric deformation matrix, ^. +A _|.[6e], and an anti -symmetric spin 
matrix, as in (5.28-5.30). 

5.3 Strain-Displacement Relations 

Strain measures valid for arbitrarily large deformations are pre- 
sented in this section. Correspond ing expressions for the variations of 
these measures are derived for subsequent use in virtual work equa- 
tions. The strain measures and variations are written in terms of the 
configuration at times 0 and t correspond ing to the T. L. and U® L. 
approaches. 

The measure of finite deformation is taken as the difference in the 
squared lengths before and after deformation of a line segment (ds) in 
the infinitesimal neighborhood of a point. 

deformation = "^ds 2 - °ds 2 = + {dx} T ^{dx} - °{dx} T °{dx} (5.33) 

The deformation may be written in terms of the initial line segment com- 
ponents, °{dx}, using the deformation Jacobian to yield 

+ ds 2 - °ds 2 = 2 °{dx} T + [el °{dx} (5.34) 

o 

t 

The symmetric. Green strain matrix, q [ e ] , appearing above is defined by 

+ [e ] = - ( + [J] T + Ui - 111) (5.35) 

o 2 o o 

Alternatively, the strain may be referred to the final line segment 
components, {dx}, using the inverse deformation Jacobian in the form 

t 2 o 2 _ t. ,T t t f 

ds - ds =2 {dx} ^.lel {dx} 


(5.36) 
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The symmetric At mans i strain matrix, ^[e], appearing above is defined by 
*[el = j ( [I 1 - °[J1 T °[J] ) (5.37) 

in which 

?[ J 3 = + [J1 _1 (5.38) 

t o 

Finally, it is verified by direct substitution that the Green and 
Almansi strains are related by 

f [e! = + [J] T *(e] + (J] (5.39) 

o o t o 

The Almansi strain matrix is often referred to as the "true" finite 
strain measure since the deformed configuration provides the reference 
state. When ^[J] represents a rigid rotation, the transformation in 
(5.39) is recognized as a simple change of reference axes. 

The symmetric strain matrix may be written in 6x1 vector form 
as ^{e} and ^{e} for the Green and Almansi definitions respectively. 
Shear strains from the matrix definition are doubled in the vector 
representation. The Almansi to Green strain transformation in (5.39) is 
written in vector form as 

%) = !tTl (5.40) 

o o t 

in which the terms of the transformat ion matrix t [T] are given in the 
appendix to this chapter (Section 5.9). 

Both the Green and Almansi strain components correctly describe the 
finite deformation due to a non homogeneous displacement field within the 
infinitesimal neighborhood of a point °P that displaces to ^"P during the 
motion. The Green strain components are employed in both the T. L. and 



U. 1. formulations. The incremental form of the Almansi strain is used 
to develop the material constitutive models. 

It is useful to separate the Almansi strain into linear and non- 
linear components by first writing 


°m » mi - In 


(5.41) 


in which 


- A fJ > 


(5.42) 


Substitution of (5.41) into (5.37) yields 


t It t T t T t 

+ [el =2 ( t lJ1 + t tj! ~ t tJI t tJ1 } 


(5.43) 


The linear displacement gradient terms are now obvious and are denoted 

by 


*[el = Ij] + (jl T ] 


(5.44) 


t 

The matrix ,[e] is often termed the linear "true" strain. The 
t 

similarity of the above expression with (5.29) is noted. If a virtual 

Almansi strain is derived from a virtual displacement field, the linear 
i* + 

term of + [5el is simply + [5e), the virtual deformation matrix. 

In the U. L. formulation it is necessary to define the strain at 
time t+At with respect to the configuration at time t. By the same 
procedure used to develop (5.33), the deformation increment is given as 


t+At , 2 
ds 


Substitution 


t 


of 


ds 


2 t+At f .T t+At f . 
= {dx} {dx} 

(5.8) for ++A+ {dx} in 


+ T t 
(dx> T {dx} 


(5.45) 


terms of the incremental 


deformation Jacobian yields 
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t+At , 2 
ds 


+ , 2 0 fr, -.T t+At, , t r , i 

- ds = 2 {dx} ^[el {dx} 


(5.46) 


in which 


t+At, , 1 ,t+At, ,,T t+At,,, fl1 , 

.[el = ( .( J 1 . [ J 1 — [II) 

i 2 t t 


(5.47) 


represents the symmetric Green strain matrix during the motion from t 
to t+At but referred to the configuration at time t. The 6x1 vector 

t+Af 

form of this strain is denoted ^{e} . This expression is exact; no 
restrictions are placed on the magnitude of the displacement increment 
^"{Au} . 

The strain defined in (5.47) cannot be added to that defined in 
(5.35) to determine the Green strain at t+At referred to time 0. 
However, using (5.8) it is easi ly shown that 


t+At r . t , t fl1 T t+At t 
[el = [el + [ J 1 . [el t J 1 

o o o to 


(5.48) 


or in vector form 


t+At, , t. . t t+At. , 
0 {s) = 0 {e} + 0 (T! t U> 


(5.49) 


This transformation implies that, relative to the Initial configuration, 
the strain increment employed in U. L. is in actuality an A I mans i strain 
increment referred to the configuration at t. 

When the body is subjected to a virtual displacement field, a 
correspond ing virtual strain field is produced. Application of the 
first variational operator, 5, on the Green strain matrix of (5.35) 
yields 

t 1 t T t t T t 

[<5e 1 = ^ ( [J1 [«JJ + [ 5 J 3 [J] ) 

o Zoo o o 


(5.50) 
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Substitution of the variation of the deformation Jacobian from (5.23) 


permits the 

above 

equation to be 

wr i tten 

i n the 

form 



f [6e ] 

o 

( 

f T t 1" 

J J] I [ 5 v] T [J] + T [J 

o t o o 

T t 

1 5v] 

Tt Ul , 
o 


(5.51) 

which may be simpl 

1 if ied using (5 

.28) to 

provide 





- + u: 

o 

T t t 

1 ! ( 5 e] [J] 

t 0 





(5.52) 

The Green 

stra i n 

rate matrix 

measured with 

respect 

to the initial 

configuration, is 

symmetric and 

is given by the transformed. 

s ymmetr i c 

deformation 

matr i 

x *[<Sel (which 

i s the 

1 i near 

portion 

of the virtual 

A! mans i strain). 

t 

Since [$el vanishes 

under a 

rigid virtual 

d ispl ace- 


ment, the virtual Green strain also vanishes. Using vector forms of the 
virtual quantities, the above transformation is written as 

*{<$£> = *[T! *{6e} (5.53) 

Inverting the above relationship yields an expression for the 
deformation vector 

!{$e} = + fT]~ 1 + {Se} (5.54) 

t o o 

Using a differential operator matrix which acts on the virtual 
displacement vector, the virtual Green strain vector is easily evaluated 
as 

+ {6e} = + [B ] + {<$u} (5.55) 

o o 

The terms of this operator matrix are defined in the appendix to this 


chapter 
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For the U. L. formulation, the virtual strain at time t+At referred 
to the configuration at time t is required. Application of the varia- 
tional operator on (5.47) followed by a substitution for the variation 
of the incremental deformation Jacabian from (5.32) yields 


t+At... , t+At. m T t+At.. , t+At,. 
t (Sel - t tJl mt [Sel f IJ 


(5.56) 


Use of the vector form of the deformation and strain matrices permits 
the above transformat ion to be written as 


t+At, \ t+At t+At. . 
t {6e} = t (T1 t+At {<Se} 


-J* 

in which |[T] is defined In the appendix. 
(5.55) for U. L. is described in Chapter 6. 


(5.57) 


An expression analagous to 


5.4 Stress Measures and Rates 


The most common stress measures employed In finite deformation 

theory and their correspond ing rates are described in this section. In 

the deformed configuration, the "true" or Cauchy stress components 

provide the natural measure. As in classical elasticity theory, these 

components are defined by considering the equilibrium of a differential 

tetrahedron extracted from the body at any time t and acted upon by the 

differential force vector ”*"{dF}. The tetrahedron has three surfaces 

parallel to the fixed reference axes. The resulting stress matrix, 
t 

denoted la\, is symmetric In the absence of body couples due to the 
orthogonality of the reference planes. 

Alternatively, stresses acting on the undeformed areas °{dA} may 

■j- 

be defined as those resisting the same differential force vector {dF}. 
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Using (5.4) , the stresses can be transformed from deformed to undeformed 
areas such that of the 

+ {dF} = f [ ct 1 + {dA} = f |j| + [ a 1 + [J]~ T °{dA} (5,58) 

in which ^[a] is the 1st Piola-Kirchoff or Lagrange stress. This matrix 
is in general nonsymmetric since the reference areas on which the 

stresses act are not orthogonal. A symmetric stress matrix is obtained 
If the actual force vector, ^"{dF}, is considered to be defined by a 

transformed pseudo force vector °{dF}. Force °{dF} is assigned a trans- 

1° f 

formation to { dF} in the same manner that {dx} transforms to 

°{dx}, i.e., 

+ {dF> = + [J] °{dF} (5.59) 

o 

Substitution of the transformed force vector into (5.58) yields 

+ {dF> = + tS ] + {dA} = f I J 1 t [J]" 1 + Cal t [J]” T f {dA} (5.60) 

in which ^[S] represents the 2nd Piola-Kirchoff (2nd P-K) or ’'pseudo" 
stress referred to areas at t=0. The above transformation shows that 
this stress matrix is symmetric; however, it cannot be physically inter- 
preted as stress in the usual manner since it resists a transformed 

force vector and does not act on deformed areas. Defining 6x1 stress 

f + 

vectors, { 0 } and ^{S}, and noting the similarity of the above 

equations with (5.39), the transformation in (5.60) may be written in 

vector form as 

I {s} = Ini l [j] ~ T +{a} (5 * 61) 

o o 1 1 o 

+ {a} = "!jjf 1 T1 T ^{S} 
o 1 1 o o 


(5.62) 
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•f 

in which matrix [ T ] was introduced In (5.39). 
o 

Now consider the configuration at time t+At in which Cauchy 

4 - 

stresses {a} act. These stresses may be referred to correspond i ng 

areas at t=0 or alternatively at time t as follows 


t+ ^ {S J 

- ++A >I 

t+At [jj-T 
0 

t+4t {.) 

(5.63) 

++4 J{S> 

- ++a Jni 

t+ >r T 

++4+ {o) 

(5.64) 


in which ^[T] was introduced in (5.57). 

Stress rates and increments are employed in the development of path 
dependent, Incrementa.l constitutive models. The 2nd P-K rate is 
obtained by differentiating (5.60) with respect to time (or the pseudo 
load parameter). The formal stress increment is obtained by multiplying 
the rate by the time increment. The differentiation process parallels 
that used in (5.505 for the variational operator. Denoting a rate (as 
distinguished from a variation) by a dot superscipt, and using (5.26- 
5.30) with a dot replacing the variational operator, the 2nd P-K stress 
rate Is given, after some manipulation, by 

f [Sl = + lJ| + |J]" 1 llo T 1 t |J]" T (5.65) 

o o 1 1 o t T o 

in which 

J[a T l = f [a 1 - + [a 1 [w] T - (wl t [<j] (5.66) 

- + tal (e] T - [el + ter] + Tr(M) + Eal 

In (5.66), [wl and [e] are at time t. The term ^ [ a ) is the actual 

t • 

Cauchy stress rate. The combined terms in (5.66), denoted ^[a^l, are 
referred to as the Truesdel I stress rate in the I iterature. Because the 
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2nd P-K stress and, therefore, its rate are symmetric, the Truesdell 

rate is symmetric by virtue of the symmetric transformation in (5.65). 

In addition, the Truesdell stress rate transforms to an equivalent 2nd 

P-K rate in the same way that the Cauchy stress transforms to 2nd P-K 

1 ” • 

stress, see (5.60). Under a rigid rotation, [SI must vanish since no 

change of the basis unit vectors occurs. The 2nd P-K stresses are thus 

invariant under a rigid rotation. This requirement necessarily implies 
t* • 

that must vanish under a rigid rotation, i.e., 

[0] = f la I - + t a] (wl T - [wl + [a! (5.67) 

Jaumann [5.7] adopted the above terms within the Truesdell stress rate 
as a suitable spin-invariant stress rate for use in constitutive laws at 
finite deformation, i.e., 

Jlffjl * + ta! - + tor I tw] T - EwI + M C5.68) 

In the case of a pure deformation increment, [w] = 0 and the Jaumann 
rate is the actual Cauchy rate. The Jaumann stress rate is also 

symmetric. In vector form, the above expressions may be written 

f ® f . . + — T + o 

J S} = l \ J r lTl (5.69) 

f • t * + t © 

f {a T } = t (aj} - + [Q] f {e} (5.70) 

1 “ 

Matrix [Q I is defined in the appendix to this chapter. The Jaumann 

stress rate finds use in finite strain plasticity as the quantity that 

is linearly related to the rate of the deformation matrix. As discussed 

by Hill [5.8], this is equivalent to a true stress vs. logarithimic 
strain relationship in simple tension. Use of the Jaumann stress rate 
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as the resultant of the constitutive model yields a zero Treusdel I and 
2nd P-K rate under rigid rotation. 

In the course of finite element computations, the equilibrium 
configuration at time t+At is sought using the equilibrium configuration 
at time t as the initial condition. Increments of the computed stress 
components are accumulated as the iterative solution converges toward 
the correct equilibrium configuration at t+At. In the T. L. approach, 
this process may be written as 

t+A+ {S} = + {S> + S {AS} (5.71) 

o o o 

As prevousiy noted, formal stress increments, {AS}, are obtained from 
the correspond ing stress rate multiplied by the time increment. At. In 
actual computations, the stress increments are determined from strain 
increments which are Induced by a displacement increment {Au}. The 
displacement increment can be thought of as a velocity multiplied by the 
time increment. 

The total Cauchy stress at time t+At may be computed using the 

transformation in (5.62). The incremental decomposition of stress above 

implies that material constitutive models produce increments of 2nd P~K 

stress. However, 2nd P-K stresses are generally unsuitable for material 

constitutive relationships which are naturally cast in terms of true 

stress- logarithmic strain. Use of a T. L. formulation necessitates some 

additional stra i n-stress transformations to obtain the {AS} needed in 

O 

(5.71). These transformations are described in Section 5.6. 

Equation (5.71) also serves as the starting point in developing an 

appropriate incremental stress decomposition for use with the U. L. 

t -1 t T 

approach. Both sides of (5.71) are multiplied by |j| [TI to yield 
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++A |{S} = + {a} + £ + {AS} (5.72) 

in which S} are 2nd P-K stresses at time t+At referred to the con- 
figuration at time t. The Truesde! i stress increments ( AS } ) appear 
naturally in this transformation. The inverse of (5.64) is applied to 
obtain the Cauchy stresses at time t+At from (5.72). The above rela- 
tionship is appealing from a physical viewpoint if the motion from t 
to t+At is a simple rigid rotation. In such a case, the accumulated 

Truesdel I increments in (5.72) are zero. The Cauchy stress 
at t+At, obtained using the inverse of (5.64), becomes transformed 

through the rigid rotation, i.e«, fTj is formed from an orthogonal 
deformation Jacobian. 

5.5 Principle of Virtual Displacements 

The principle of virtual displacements is applied to derive the 
equilibrium conditions for a body undergoing arbitrarily large displace- 
ments. The basis of the argument is that, for an imposed virtual 
displacement field, the virtual work remains invariant with respect to 
the configuration in which the variables are measured [5.9]. 

Consider a body at time t+At which occupies a volume and is 

bounded by the surface area ++A+ A. Let the body be in equilibrium under 
a set of body forces { p} and surface tractions {q}. For dynamic 
analysis, the inertia effects are Included in the body forces according 
to D'Alembert*s principle. If a virtual displacement field {<Su} 
consistent with the kinematic boundary conditions is applied to the 
body, the virtual work done is given by (ignoring thermodynamic effects) 



i 90 


t+At 


SW 


/ 

t+At 


t+At 


C5u> T t+4+ {pi t+4+ dV 


t+At 


,T t+At, t+At,. 
{5u} ++A+ {q> dA 


A 


(5.73) 


and must be zero. The mass density is denoted by p. The body forces 
are in terms of unit mass; the surface tractions are in terms of the 
force per unit area at time t+At. The applied surface tractions {q} are 
written in terms of the equilibrating Cauchy stresses. The classical 
procedure [5.9] is then employed to convert these integrals to the 
virtual work form using Gauss* theorem. The result is 


++A+ SW = ++A+ SW + ++A+ SW. = 0 

ext i nt 


(5.74) 


in which the external virtual work is given by 


t+4t, w _ r t+At , ,T t+At . . t+At 

5W axt ’ W„ 4 tSuJ o lp) dV 


(5.75) 


+ U, tSu)T dA 


The internal virtual work is given by 


t+At 


SW, 


int 


t+At. 


Tr( 


t+At 


la] 


t+At, 


Sel ) 


t+At 


dV 


(5.76) 


The internal virtual work measured from the instantaneous configuration 
involves the Cauchy stress and the variation of the deformation matrix 
(which is the linear part of the Almansi strain variation). The above 
expression may be simplified using vector forms of the strain and stress 
to yield 


t+At 


SW. 


int 


= " I 


t+At 


t+At, 


{5e} 


T 


t+At, ,t+At ... 
{ 0 } dV 


(5.77) 


To apply this equation, the terms must be expressed as functions of 
variables in a known equilibrium configuration. The obvious choices are 



the configuration at time 0 (T. L.) and time t (U. L. ) . To develop 
virtual work expressions for the Total Lagrangian formulation, all 
variables are referred to the undeformed configuration. Use is made of 
the following transformations previously derived 


tt&t 

t+At 

t+At 


{a} 
{Se} 
dV = 


t+ 4 t m T t+ 4 + (s) 

o 1 1 o o 


ttAtr-.-l t+At,. , 
[T] {6e} 

o o 


t+A >l °dV 


(5.78) 

(5.79) 

(5.80) 


Upon direct substitution of these transformations into (5.76), the 
internal virtual work becomes 


t+At 


<5W. 


nt 



m + T wt o dv 

O O 


(5.81) 


The internal virtual work is In terms of the 2nd Piol a-K i rchof f stresses 
and the variation of the Green strain. In (5.9] the above internal 
virtual work expression is derived by direct transformation of the 
integrals in (5.76) to time 0 before Gauss ! s theorem is applied. 

To obtain the U. L. formulation, the terms in (5.76) are written as 
functions of the known equilibrium configuration at time t. The pre- 
viously derived transformations are 


t+At. . t+At. ,,-1 t+At T t+At j. , 

{a} = + |J| + [T] + {S} 

(5.82) 

++4t Ue} • t+4 *rrr’ 

t+At.. . 
+ {se} 

(5.83) 

++4+ dV - t+4 ||j| f dV 


(5.84) 

Upon substitution of these 

transformations 

into (5.76), the internal 


virtual work expression becomes 
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++A+ SW. , = - /. ++A t(5eT ++A t{S> + dV 
t w + t 


i nt 


V 


+ 


(5.85) 


which is analogous to (5.81) except that integration is performed over 
the configuration at time t rather than time 0. 

Transformation of the external virtual work integrals (5.75) to 
either T. L. or U. L. forms introduces another level of complexity. How- 
ever, the volume integral for the body forces presents no difficulties 
and is readily transformed to 


t+At ,BF 


SW° r = / °p {Su} 

ext J o„ 


T t+At, . o 

Ip} dV 


(5.86) 


for T. L. and to 


t+At... BF f t ,T t+At. . t ... 
5W ext = J p {p} dV 


(5.87) 


for U. L. These transformations are particularly simple since the di- 
rections of the body force components are not a.I tered under the differ- 
ential volume conversions. Body forces are conservative by defini- 
tion. Inertia effects are simple to incorporate in (5.86) by the sub- 
stitution of {-u} for {p}. Similarly, for U. L. , inertia effects are 
i ncorporated by transforming the integral in (5.87) to time 0 by noting 
that ^p^dV = °p°dV (conservation of mass) and letting {p} = {— u> . The 
transformed integral is Identical to (5.86). This is of major practical 
consequence as the mass matrix for the U. L. formulation is identical to 
that for T. L, and remains constant for all times throughout the 
response . 

The transformation of external virtual work for the applied surface 
tractions in (5.75) is more complex due to the different orientation 
of d " +A '*"dA when referred to times 0 and t. This dependence of external 
work on the displacements leads to non-conservative loading. An 
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interesting case of applied surface tractions is a normal pressure 
loading as considered by Oden [5.101, Nayak [5.1], and Mackay [5.111. 
Thus, 


t+At, , t+At t+At t+At, .t+At 

++A+ {q> dA = q n {n} dA 


(5.88) 


where ^" 4A ^{n} contains the components of a unit outward normal to the 
deformed area ^" 4A ^"dA and ^ +A ^"q n is the intensity of the appl ied normal 
pressure (suction is positive). Transforming (5.88) to time 0 using 
(5.4) yields 


t+At t+At , , t+At,. t+At t+At | ,| t+At f ,, -T o, , o.„ on . 

q t n} dA = q J [J] {n > dA (5.89) 

n n o‘ 1 o 


where °{n} is the unit outward normal at the same material point on the 
undeformed element surface. The surface traction integral of (5.75) 
becomes 


t+4t ^ t 


f {Su} T ++At [Jj' 
°A ° 


•T o, 


t+At, 

0 |j 


dA 


(5.90) 


in which °{q*} has been substituted for ++At q n °{n>. The analogous 
expression for Updated Lagrangian is 


t+At.,, ST f .T t+At. . ,-T . t+At i , t.. 

AW ex+ = J + {6u> + [JJ f (q*} + |J| dA 


(5.91) 


« t t+At t 

with {q*} = q^ {n}» When the displacements represent large rota- 
tions and/or finite geometry changes, the use of (5.90) or (5.91) is 
essential to correctly generate residual loads in finite element 
computations. When only concentrated (non-follower) forces and body 


forces are applied, the surface Integral above vanishes. Similarly, if 
the pressure loading is applied on a surface that undergoes 
infinitesimal displacements, equation (5.91) degenerates to (5.90) in 
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mechanics type problems with loading far removed from the intensely 
deformed plastic zone near the crack tip. 


5 . 6 Constitutive Models for Finite Deformation El asto-P I ast i c ? ty 


Constitutive models provide a connection between incremental strain 
and stress changes. There exist no generally applicable constitutive 
relations for materials subjected to finite strain magnitudes. However, 
two special cases have received considerable attention by finite element 
researchers and have resulted in reasonable engineering approximations for 
material behavior. 

Oden 15.10] has reviewed formulations for nonlinear elastic (hyper- 
elastic) materials similar to rubber. For such materials, it may be 
possible to develop an elastic energy function in terms of the total Green 
strain components. The 2nd P-K stress is obtained by partial differentia- 
tion of the energy function with respect to strain. Thus, 

_ +• 

W = W ( {x}, {s}, material constants) (5.92) 
o 


>) . 

O „t 


(5.93) 


a'{e) 

o 


The relationship between increments of strain and stress is obtained by 
again differentiating (5.93) which yields 


d ^{S} = *[DI d Ue} 
o o o 


(5.94) 


in which 
t 


o 


[Dl = 


3 2 W 


3*{e} ^{e} T 
o o 


(5.95) 
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The matrix ^[D! contains the tangent elastic moduli. When the material is 
isotropic, the function W may be simplified through the use of strain in- 
variants. The commonly used Mooney-Rivl in model for rubber assumes incom- 
pressible behavior with an energy function defined by the first and second 
strain invariants and two material constants. 

The hyper-elastic material models are of little use for most en- 
gineering materials which experience permanent deformation when subjected 
to stress states outside an elastic domain. The remainder of this section 
focuses on an approximate theory for finite deformation, elasto- 
plasticity. This theory has been adopted by a number of finite element 
researchers, especially for applications involving the elastic-plastic 
fracture of metals. The fundamental assumptions of this theory are; 

1. the material is Initially isotropic; 

2. the materia! work hardens isotropically; 

3. there exists a linear relation between stress and strain 
i ncrements; 

4. additive decomposition of strain increments into elastic and 
plastic parts is valid; 

5. recoverable elastic strains are infinitesimal; 

6. the stress increment is independent of the rate of rigid body 
rotation. 

The basis for this material modeling theory is provided In Hill's 
work (5.81. More recently, Nemat-Nasser [5.12] reviewed the status of 
developments of finite deformation plasticity and provided additional 
arguments In support of assumption (4). His argument is that 
infinitesimal incremental deformations within a crystalline solid 
occurring from a finitely deformed configuration must satisfy additivity 
as in the infinitesimal deformation case. He further argued that the 

f e 

symmetric part of the velocity gradients (the deformation vector ^{e}) may 
be decomposed into elastic and plastic components. 
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Given the above restrictions, the Prandtl -Reuss equations are adopted 
as usual, but in terms of true stress and true strain. Strain rate ef- 
fects, creep, and thermal dependent material constants may be incorporated 
in this theory following the same techniques used in small strain 
plasticity, although little of this work has been reported in the 
I iterature. 

Following Hill [5.8], Nemat-Nasser [5.12] and Hibbitt, et al. [5.13], 
the Jaumann stress increment and the increment of the deformation vector 
provide appropriate stress and strain increments which vanish under rigid 
body motion. The incremental relation becomes 

^{Aaj} = J[D t 1 *{Ae} (5.96) 

in which ^[Dy] is formed exactly as in the case of infinitesimal strain 
plasticity [5.14], but using true, Cauchy stresses. The above relation- 
ship is strictly valid for differential changes of the deformation vector. 
The von Mises yield criterion and associated flow rule are adopted for 
metals. For stress states inside the yield surface, the tangential con- 
stitutive matrix, ^[Dy], is composed of simple elastic constants. For an 

■j* 

elastic-plastic state, the terms of y[D] are derived from the elastic con- 
stants, the current stress state, the material strain hardening character- 
istics, and the history of Cauchy stresses. For uniaxia! tension, it is 
simple to show that (5.96) represents a true stress vs. logarithmic strain 
relationship. Numerical refinements, such as the subincrement method, 
used in small strain plasticity are equally applicable in (5.96) to assure 
satisfaction of the flow rule. 

The incremental relation in (5.96) is adopted in both U. L. and T. L. 
formulations through appropriate transformation of the terms to the re- 
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quired reference conf igura i ton. Considering first the T. L. formulation, 
the deformation increment is obtained from the Green strain increment 
using (5.54). 

*{Ae} = Jm" 1 £{Ae} (5.97) 

The assumption here is that the Green strain increment, ^{Ae}, is actually 
of differential, not finite, magnitude (see equations 5.50 - 5.55). The 
resulting Jaumann stress increments obtained from (5.96) are transformed 
to Truesdel I increments with (5.70). 

^{A 0 t } = ^{AOj} - f f lQ] *{Ae} (5.98) 

The Truesdel I increment is transformed to the required 2nd P-K increment 
with (5.69). 

+ {AS} = + I J | f [J]~ T J{A0 t > (5.99) 

To compute a tangential stiffness, matrix, it is necessary to have a 
direct relationship between increments of 2nd P-K stress and Green strain 
in the form 

+ {AS} = + [D t ] + {Ae} (5.100) 

O O I o 

Combining (5.96-5.99), the tangent material matrix is given by 

*[D T 1 = + |J| Vf T [ + (0 T ! - *[Q1 ] Vf 1 (5.101) 

oT o 1 1 o '-TT t J o 

However, because ^[Q] is not symmetric, the resulting constitutive rela- 
tion is non-symmetr ic which leads to an undesirable non-symmetr ic stiff- 
ness matrix. Osias (5.151 arrived at this same relationship (5.101) and 
retained the non-symmetr ic terms. Upon examination of the [Q I matrix, the 
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non-symmetr ic terms are found to arise from the trace of the deformation 
matrix which represnts the incremental volume change, see Eq. (5.66). If 
these terms are neglected, then ^[Q] becomes symmetric. Recalling 
assumption (5), the rationale to neglect these terms is that: (1) under 
elastic conditions, the volume change multiplied by the stresses is a very 
small term compared to the elastic moduli, and (2) under perfect plas- 
ticity conditions the incompressibility constraint forces a zero volume 
change. McMeeking and Rice [5.16] omitted the trace term as did Nagtegaal 
and de Jong [5.17], although neither group noted this assumption. Further 
complications arise when the material has a very low hardening modulus (in 

terms of true stress- logarithmic strain). The intermediate result, 
t t 

^[D^J - [ Q ] , may develop negative terms on the diagonal for large tensile 

1 " 

strains. One possible remedy is to simply omit + [Q] in (5.101) when 
forming a suitable constitutive relation for use in computing a tangent 
stiffness. The writers have demonstrated that this technique is 
acceptable provided the exact stress increment is computed using (5.96- 
5.99). 

Consider now constitutive relations for use with the U. L. formula- 
tion. The deformation increment at time t+At for use in (5.96) is com- 
puted from the Green strain Increment referred to the configuration at 
time t using (5.57) . 


ttAt 

t+At 


{Ae} 


t+At -1 t+At 
+ m + {Ae} 


(5.102) 


Again, it is implied that + {Ae} is actually of differential magnitude. 
The Jaumann stress increment is obtained from (5.96) and transformed to a 


Truesdel I increment by 


t+At 

t+At 


{Acf.j.} 



(5.103) 
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The Truesdel I increment is transformed to an increment of 2nd P-K stress 
referred to the configuration at time t (which is just the Truesdel i 
stress transformed to time t) 

t+At r . c \ _ t +At | . | t+At r T , — T t+At,. , ,c 

+ {AS} = t |J| f f T ] ++At Ua T ) (5.104) 

In U. L., the tangent stiffness matrix is formed for the configuration at 
time t and generally held constant during the increment At. Thus, the 
material tangent matrix is given by 

Jid* t 1 = J[D t 1 - JlQl (5.105) 


in which both matrices on the right side are generated in terms of actual 

•f 

Cauchy stresses at time t. The [D*l above relates increments of Green 
strain and Truesdel I stress at time t. If the stiffness is updated 
between t and ttAt (to accelerate convergence) the new constitutive matrix 
is evaluated as in (5.105) using the most current values of Cauchy 
stress. The same difficulty that arises due to the non-symmetry of [QJ 
for T. L. also occurs here for U. L. 

Approximating assumptions for (5.96-5.105) have been introduced by 
various Investigators as summarized below. Altur?, et. al. [5.18! 
employed an U. L. approach with the approximation 


t+At, , _ t+At t+At. , 
t {AS} ~ t+At °T t {A£} 


(5.106) 


which neglects changes in the structural configuration over the inter- 
val At. In forming the constitutive matrix, [Q] was simply ignored. 
Hibbitt, Marcal, and Rice [5.13! fully developed the T. L. equations 
(5.96-5.101) but implemented a simplified form valid only for small 
strains by neglecting IQ I in a I I equations. The determinant of the 
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deformation Jacobian was also assumed to remain unity which effectively 
reduces [T] to a rotational transformation. Based on these assumptions, 
they claim that the implemented formulation is valid for the case of large 
rotations but infinitesimal strains. McMeeking and Rice [5.16] developed 
an U. L. approach in which the configuration is continuously updated. 
Under such a procedure, the transformation in (5.102) is unnecessary, 
i.e., [T] is an identity matrix. In addition, the transformation of 
Truesdel! stress in (5.104) becomes unnecessary. Bathe, et al. [5.5] 
presented both T. L. and U. L. formulations. For T. L. , the material 
constitutive relation in (5.100) was formulated using classical plasticity 
theory but in terms of total 2nd P-K stresses. The proposed U. L. 
constitutive model (5.106) is that adopted by Alturi, et a I . [5.18]. 

With exception of the Jaumann to Truesdel I transformation in (5.98), 
differences in the various U. L. approaches vanish for sufficiently small 
motion over the interval At. However, the choice of a "suff Iciently" 
small At cannot be assured for complex structures undergoing large dis- 
placements; thus, the more exact transformation In (5.104) is preferred. 
With this approach, the T. L. formulation is recovered as the limiting 

case when no updates of the deformed configuration are performed. The use 
+ 

of (5.100) with q [Dj] formed in terms of 2nd P-K stresses, as proposed in 
[5.5], is clearly unacceptable for general applications. 

There continues to exist some question as to the need for the [Q] 
terms appearing in the various transformations. These terms arise from 
the Jaumann to Truesdel I stress increment transformat ion. In the elastic 
range, these terms are on the order of stress, which when multiplied by a 
strain increment are clearly negligible compared with elastic moduli. In 
the plastic range, the terms of [Q 1 become more comparable in magnitude to 
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the plastic moduli (which may be on the same order as the stress) unless 
the deformation increment is truly of differential magnitude. Hibbitt, et 
al. [5.13] acknowledges the role of [ Q I then chose to ignore it under the 
restructions of small strain-large rotation. McMeeking and Rice [5.163 
include [Q] in their continuously updated formulation as an initial stress 
matrix and do not report any numerical difficulties. They conclude that 
whenever the slope of the Cauchy (true) stress vs. logarithmic strain 
curve has a magnitude comparable to the current stress level, the 
predicted tangent stiffness cannot be accurate. Both [5.13] and [5.161 
conclude that the relative importance of [Q] terms is not fully known and 
requires further study. 

5.7 Summary and Comparisons of the U.L. and T.L. Formulations 

Using matrix notation, the equations of nonlinear continuum mechanics 
that provide a rational basis for finite element analysis have been 
presented in this chapter. Two formulations are described in detail; 
namely, the Total Lagrangian (T. L.) and Updated Lagrangian (U. L.) The 
initial configuration of the body at time t=0 serves as the reference 
state for all variables in the T. L. formulation. The configuration at 
time t serves as the reference state to describe all variables at time 
t+At in the U» L. formulation. Full details of the stress measures, 
rates, transformations and nonlinear constitutive models have been 
described that enable analyses involving both infinitesimal and finite 
strain magnitudes. Finite magnitude strains commonly appear in structures 
that experience transient response after being subjected to localized 
blast or impact loadings. 
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The two nonlinear formulations, T. L. and U. L. , are shown to derive 
from a common definition of the rate of work per unit mass which leads to 
equivalent virtual work expressions when the correspond i ng reference 
states are introduced. Analyses derived from each formulation are iden- 
tical to within truncation errors introduced in the numerical proce- 
dures. However, as noted throughout this chapter, various investigators 
have introduced simplifying approximations into the formulations that 
limit their applicability to the most general class of problems. The U. 
1. formulation offers the most temptation to simplify the equations since 
one might argue that the difference in configuration between time step t 
arid t+At is sufficiently small that all second order nonlinear terms may 
be discarded. Stress rates and transformations in the correspond i ng 
constitutive models are simplified similarly. No approximations of any 
type are permissible in the T. L. formulation since the deformed and 
reference configuration are seldom similar. 

The differences in computational efficiency of the two formulations 
lies in the stress rate transformations and the strain-displacement rela- 
tions. The mass representation is identical in both formulations and thus 
does not enter into the discussion. Similarly, computations associated 
with a nonlinear constitutive model are identical as the computed results 
required by the model are identical in both formulations. As shown in 
Section 5.6, the same number and type of transformations to account for 
stress rates are required In each formulation. The U. L. approach appears 
to gain a computational advantage if the simplifying approximation is 
employed that eliminates the Truesdel I to 2nd P-K stress transformation. 
However, this advantage is offset by two factors: (1) the number of time 
steps to assure that the change in configuration over a single time step 
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is "small" remains unknown, and (2) the stresses at the beginning of the 
step and the increments that occur during the step must be transformed to 
the configuration at the end of the step. This last process is unneces- 
sary for the T. L. formulation in which all stress increments are trans- 
formed to the initial configuration as computed. The stress increments 
are simply added to the existing stress state. Therefore, in terms of 
stress computation, the T. L. approach actually has a slight efficiency 
advantage over the exact U. L. formulation. 

The computation of an increment of the deformation vector is accom- 
plished in the same symbolic form (5.53 and 5.57) for both formulations. 
The differences lie in the effort required to construct the [T I matrix in 
each case. For U. L., the formulation of derivatives in £ T] requires that 
updated coordinates be available, whereas in T. L., derivatives are always 
formed with respect to the initial configuration. These derivatives could 
be computed once and saved for re-use. Some effort is necessary to 
continuously update the coordinates. More importantly, the strain incre- 
ments produced by U. L. cannot be simply accumulated to form the total 
strain, as shown in (5.48). This transformation to a common reference 
configuration requires additional computational effort that is unnecessary 
in the T. L. formulation. 

In summary, an examination of the computational efficiency of the 
transformations necessary in each approach, reveals that T. L. appears to 
have a slight advantage over U. L. This does not include any factors that 
are Introduced when the governing equations are cast into a form suitable 
for finite element analysis. In a practical sense, the advantages of 
either the U. L. or T. L. continuum theory appear computationally 
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insignificant given other costs of computation, for example, equation 
solving and massive data transfers between memory and disk. 
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5. 9 Appendix — Notation 


The following convention for vector and matrix subscripts and super- 
scripts is used: 

1. a left superscript indicates the discrete time configuration In 
which the variable occurs, 

2. a left subscript in conjunction with a left superscript indi- 

cates the discrete time configuration with respect to which the 
variable is measured, 

3. a left subscript by itself indicates an Increment from time t 

to t+At referred to 'the configuration at the specified time, 

4. a dot over a symbol denotes a rate quantity. Left subscripts 
and superscripts indicate the time and reference configuration, 

5. a "5” symbol denotes a variation of the quantity. Left sub- 
scripts and superscripts indicate the time and reference 

configuration. 


[ I 


I ! T , { } T 

°{x), + {x}, ++A+ {x} 


{ !> 

f [Jl 

o 

tj. , ttA+j. , 
{u}, {u} 


Square brackets denote a matrix 
Curly braces denote a column vector 
Right superscript "T" denotes the transpose 
Cartesian coordinates at time 0, t, and t+At 
unit vectors 

Deformation Jacobian matrix 

Vectors of displacements from 0 to t and from 0 



o t 
P. P 



to t+At 

Vector components of the surface normal at times 
0, t 

Mass density at times 0, t 
Area at times 0, t 
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*(q>, 

o {SJ - o ,SI 


t+At. . t+At 
t 1 ^' t 


• "t’ • 

f £cr T > , f [ a T 1 


t , • , + f • , 

+ l«jl 


t+At f . t+At 

1 £ / p 


t+At r , t+At 
„{e}, 

o o 


[ V ] 


.lei 


Jw3 


[ T 1 


[SI 


[el 

tel 


Volume at time 0 and t 

Vectors of surface force components at time t per 
unit area at 0 andd t 

Vector of body force components per unit mass at 

time t referenced to configuration at 0 

Cauchy stress vector and matrix at time t+At. 

(Note that (crl = 4x ,.[a]). 

t+At 

2nd Piola-Kirchoff stress vector and matrix at 

time t referenced to the configuration at 0 

2nd Piola-Kirchoff stress vector and matrix at 

time t+At, referred to configuration at t 

Truesdel I stress vector and matrix at time t 

referred to the configuration at t 

Rates of Jaumann stress vector and matrix at time 

t referred to the configuration at t 

A I man si strain vector and matrix at time t+At 

(note that ++A+ [ e ] = j+^tel). 

Green strain vector and matrix at time t+At 
referenced to configuration at 0 
Incremental displacement gradient matrix at t 
referenced to configuration t 

Deformation matrix at t referred to configuration 
t (symmetric part of Jjvl) 

Spin matrix at t referred to conf iguration t 
(anti-symmetr ic part of ^[v]) 

Transformation matrix to convert Almansi strains 


at time t to Green strains 
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Tr( ) 
^Q1 


<$W 

W 


t 

t 


[D t I 


Displacement Jacobian at time t referred to 

configuration at 0 

Differential operator to yield Green strain 
increments from displacement increments 
Denotes the trace of a square matrix (sum of 
diagonal terms) 

Transformat ion matrix that converts Jaumann stress 

increments at time t to Truesdel I stress 

increments at time t 

Virtual work quai ity 

Elastic energy density function 

Elasticity matrix at time t referred to 

configuration at 0 

Elastic-plastic matrix at time t referred to 
configuration at t 



2! 0 




2 1 
J n ; 

J 2 

J 21 

i 

+ m = 

o 

j 2 1 

J 12 i 

J 2 

J 22 

1 J J 
| J 1 2 J 22 


_ 2J 1 1 J 12 i 

2J 21 J 22 

' J 11 J 12 + J 12 J 21_ 



6 + {e} = [B36{u> 

o o 
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CHAPTER 6 


NONLINEAR FINITE ELEMENT EQUATIONS 


6« 1 General 

The finite element concept of discretization is combined with the 
continuum mechanics theory in matrix form in this chapter. With the 
addition of solution procedures for nonlinear transient equations, a 
very general analysis capability results. The form of the finite 
element equations of motion are first expressed without regard to a 
particular choice of reference configuration. The use of such a general 
form makes the detailed discussion of nonlinear transient solution 
procedures equally applicable for an arbitrary reference 
configuration. Specific forms of the element' stiffness matrices and 
internal resisting force vectors are derived for the Total and Updated 
Lagrangian approaches. The relative computational effort required for 
each approach and the implications of a substructured modeling procedure 
are also examined. Recommendations are made for appropriate 
formulations and solution procedures necessary to support the analysis 
of a broad class of problems. 

A number of investigators have contributed to the current state of 
rigorous finite element formulations derived from the continuum theory 
of the previous chapter. Three approaches appeared almost 
simultaneously in the I iterature of the early 1970 ! s. Hibbitt, Marcal, 
and Rice [6.1] first presented a comprehensive formulation using the 
configuration at time t=0 as the reference state (Total Lagrangian). 
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Matrices required for finite element analysis, including the non- 

symmetric stiffness for displacement dependent loading, were 
described. Using matrix notation throughout, Nayak [6.2] independently 
derived the same element matrices. 3oth of these studies addressed only 
static analysis and employed the symmetric 2nd P-K stress and Green 
strain as conjugate measures in the configuration at time t=0. Nemat- 
Nasser and Shatoff (6.31 used a combination of current (time t) and 
initial reference configurations with the non-symmetr ic 1st P-K stress 
measure. 

Investigators advocating the use of the current configuration 
(convected coordinates) for the reference state included Hartzman and 
Hutchinson [6.4], Belytschko and Hsieh [6.51, and Key 16.61. In each 

of these studies, the motion of the body over At was decomposed into a 
rigid motion and a pure deformation. The procedures worked well in each 
case even though different schemes were devised for the decomposition 
procedure. Only Key [6.61 focused on the finite strain case; the other 
two studies were concerned primarily with targe rotation effects. 

Interestingly, each of these studies addressed the problem of nonlinear 
wave propagation using explicit integration procedures. McMeeking and 
Rice (6.71 presented a rigorous formulation based on a continuously 
updated reference configuration applicable for finite deformation and 
rotat i on . 

Yaghmai and Popov [6.81 were the first to describe an Updated 

Lagrangian approach which attempts to combine the Total Lagrangian and 
the updated coordinate approaches. This was done by adopting the 
configuration at any time t prior to ttAt as a reference state. The 
existing Green strains and 2nd P-K stresses in the reference 



configuration are simply treated as initial strains and stresses to 
which are added increments of 2nd P-K stress and Green strain. This 
approach captures the spirit of a convected coordinate approach but is 
more appealing for four reasons: (1) the reference conf igurat ion does 

not need to be continuously updated, (2) the decomposition of motion 
into deformation and rotation follows Lagrangian mechanics, (3) the 
Total Lagrangian approach is recovered exactly if the reference 
configuration is not updated, and (4) no limits on the deformation 
magnitude over At are implied due to linearization if the exact U. L. 
equations of the previous chapter are utilized. Nagarajan and Popov 
[6.91 subsequently used the U. L. approach to study viscoplastic 
response of thin shell structures. Bathe, et. al. [6.10, 6.111 
formalized the U. L. procedure for general transient analysis and 
published a large number of papers using the procedure. 


6.2 Nonlinear Equations of Motion 

The virtual work principles derived in Section 5.4 provide the basis 
to generate approximate equations of motion using finite element con- 
cepts. Considering the motion of the body over the interval from t to 
t+At, and without regard to a particular reference coordinate system, 
the principal of virtual work provides that at t+At (see 5.74-5.76) 

SW = <SW , + 6W„ , = 0 (6.1) 

ext mt 

where the t+At left superscript has been dropped in this section for 
brevity. All quantities are at time t+At unless otherwise noted. 
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The internal and external virtual work terms are given symbolically 
by 

<5W. nt = -/ y *{<Se} T *{a} dV (6.2) 

6W ext = / v *p{«Su} T {p}*dV + /^ A (5u} T ^{q>*dA (6.3) 

The aster ik (*) is adopted to indicate an unspecified reference 

configuration for the variables. As shown in Section (5.4), the general 

forms of the virtual work equations for both T. L. and U. L. are 

identical. The differences derive from the choice of reference system 

(*) and the requirement of conjugate strain and stress definitions. 

This notation permits a completely general discussion of the finite 

element process and solution techniques prior to introducing specific 

matrices for the T. L. and U. L. approaches. 

The finite element concept is invoked at this point to provide a 

spatial discretization of the structure. Individual elements are 

connected at discrete points termed nodes. At time t+At, the continuous 

displacement field, {u}, within each element is approximated by a set 

of interpol ation functions, [ N I , which act on the nodal 

0 

displacements {a }. Thus, 


{u} = INI {a 9 } 


(6.4) 


In 2-D problems, {u} is a 2 x 1 vector containing u and v displacement 

components as a function of position within the element. Vec- 
0 

tor {a } then contains (2 * number of element nodes) terms. Terms of 
[ N 3 are usually simple functions of position within the element and do 
not depend on its deformed shape or time. 
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To facilitate the development of elements with initially curved 
edges and faces, the shape functions are expressed in terms of an 
intrinsic, convected coordinate system. A one-to-one mapping from the 
parent configuration (usually square) to the actual shape in the 
structure is provided by a set of interpolating functions and the 
Cartesian coordinates of the element nodes. As will be shown, the use 
of initial or updated Cartesian coordinates for this mapping plays a 
major role in the differences between T. L. and U. L. When 

interpolation functions for the element mapping and those for the 

displacement interpolat ion are identical (6.4), the popular 
i soparametr ic element family is obtained. Fortunately, the details of 
transformations between the Intrinsic isoparametric coordinate and 
global Cartesian coordinates do not affect the basic finite element 
solution procedure. 

At time t+At virtual displacements, {<Sa }, consistent with the 
kinematic constraints are imposed on the nodes. Correspond ing 
variations of the element displacement field, {Su}, are given by 

{<5u} = IN] {6a®} (6.5) 

Virtual strains caused by the virtual displacement field are denoted 
*{6e}, in which the definition of strain appropriate for the selected 
reference configuration is implied. However, it is completely general 
to express the virtual strains as 


*{5e} = L! {6u} 


( 6 . 6 ) 
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in which *[!_] is a differential operator matrix that may be a nonlinear 
function of the deformed element configuration at time t+At. Combining 
(6.5) and (6.6) yields 

*{5e} = *IL1 [N1 {<5a e } = *(B] {<5a e } (6.7) 

in which the conventional notation *(B1 of the virtual strain-nodal dis- 
placement relationship is introduced. The general form of the internal 
virtual work for an element is written 

SW e . n+ = -{Sa e } T /^ v *[Bl T *{cr}*dV e = -{5a e } T {F®} (6.8) 

in which the notation {F°} is introduced to denote equivalent nodal 
forces induced by the deformed element. The (*) left subscript is not 
required on such terms since these forces are directed along global 
coordinate axes. 

Consider now the virtual work of body forces within an element, 

«Wg F = / #v *p {6u} T {p}*dV 0 (6.9) 

in which *p Is the mass density and {p} are the Cartesian components of 
body force per unit mass. Two types of body forces are considered: (1) 
self-weight and centrifugal, In which the body force intensity is 
written in terms of current density as *{F}; (2) inertia effects using 
D’Alembert's principle with {u} defining the element acceleration 
field. Using (6.4), the body force virtual work may be written in the 


form 


6W® F = {5a e } T j^ v tN] T *{F}*dV e - {5a e } T /* v *p[Nl T {u}*dV 0 


( 6 . 10 ) 
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To simplify this expression, let {BF e } represent the first integral 
(equivalent nodal loads due to body force). Again the (*) left 
subscript can be omitted. The second integral is simplified by noting 
that in analogy with (6.4) 

{u} = IN] {a 6 } (6.11) 

•and that the second integral of (6.10) must be invariant with the 
deformed configuration (conservation of mass— see Section 5.4). The 
virtual work of element inertia forces may thus be written in the form 

6W 8 = -{6a e } T tM e ] {a 8 } (6.12) 

In which 

[M® 1 = / v [N] T [N]°p°dV (6.13) 

o v 

0 

Matrix [M ] is usually termed the consistent element mass; it is 
computed once at the beginning of the solution for the initial 
conf iguration and recalled from secondary storage whenever needed. 

In a similar manner, the virtual work of external tractions applied 
over the element surfaces may be written in the form 

<SW® T - {Sa e } T /^ A INl T *{q}*dA « {SaVcTF 8 } (6.14) 

The total virtual work for an element is thus written 

<$W 8 = {<Sa 8 } T j-{F 8 } + {BF 8 } + {TF 8 } - EM 0 ]{a @ }} (6.15) 

Virtual work for the complete structure is obtained by summing the 
contributions from each element. Element nodal displacements are 
related to the structural nodal displacement vector through a simple 



220 


Boolean connectivity matrix, which accomplishes the symbolic assembly 
process. The resulting virtual work expression is analogous to (6.15) 
except that vectors are of structural, rather than elemental, size. 

<SW = 0 = {6a} T j-{F} + {BF} + {TF} - (6.16) 

The choice of nodal virtual displacements is arbitrary; thus the summed 
terms In { } must vanish for dynamic equilibrium at time t+At. The 
result is a set of nonlinear, simultaneous equations in the nodal 
displacements and accelerations which may be written in the form 

{R} - {P} - {F} - [Ml {a} = {0} • (6.17) 

Vectors {F} and {P} are implicit functions of the nodal displacements 
and generally cannot be written as a matrix multiplied into the current 
displacement vector. For simplicity, the external applied load effects 
are combined into {P}. Likewise, damping effects have been neglected; 
however they could be included by a term of the form [Cl {a} analagous to 
the inertia term. 

At this point it is instructive to review the sources of 
nonlinearity in the dynamic equations of motion (6.17): 

1) the change in surface orientation and magnitude under loading 
causes {P} to become a function of the displacements; 

2) the Internal resisting force vector, {F}, may be nonlinear due 
to the material constitutive relationship between total stress 
and strain (*{a} and *{e}); 

3) the internal resisting force vector may be nonlinear due to the 
dependence of *[B] on the nodal displacements. 

The following section considers general solution procedures. 
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6.3 



At any given time, 0, At, 2At,.....t, t+At the nodal displace- 
ments, velocities, and accelerations are sought such that the nonlinear 
equilibrium equations represented by (6.17) are satisfied. A temporal 
integration operator is employed for the time discretization. Both 
implicit and explicit operators are applicable. With explicit 
operators, computations are performed directly on the vector form of the 
dynamic equations in (6.17) with {R} taken as zero. Component terms of 
each vector are computed directly from the integrals in (6.8, 6.10 and 
6.12). No stiffness matrices or mass matrices are ever assembled. The 
advantage of an explicit procedure is that it provides the computational 
efficiency needed to render feasible the solution of wave propagation 
problems (which require exceptionally small At). The computations over 
each time step are relatively simple. The major disadvantage is that 
nonlinear effects due to the spatial variables cannot be "iterated" out 
at constant dynamic load, i.e., no equilibrium iterations within a time 
step are possible. 

With implicit integration operators, the displacement increment over 
At is used to predict the accelerations at ttAt. The displacement 
increment is obtained by using an effective stiffness matrix to form a 
set of simultaneous equations. The advantage of an implicit approach 
lies in the capability to completely correct for nonlinear effects in 
the spatial variables that occur over At through equilibrium 
iterations. In addition, corrections in the estimated displacement 
increment over At during Iterations are utilized to improve the 
predicted acceleration at t+At. At t+At, the dynamic equilibrium 
equations can be satisfied to within a specified tolerance with an 
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implicit procedure. The major disadvantage of the implicit scheme is 
the increased computational effort required to assemble and triangulate 
the effective stiffness matrix. Thus, a solution procedure based on an 
implicit scheme is more appropriate for problems in which wave 
propogation is not as important as gross inertia effects. Much larger 
time steps are permitted thus reducing computational costs. If 
necessary, the implicit scheme can be used to compute wave propagation 
effects. 

The implicit solution scheme is adopted for detailed discussion here 
and will be adopted in the software. The implicit scheme builds upon 
the experience gained in the solution of static nonlinear problems. 
Static analysis is recovered as the degenerate case when the mass matrix 
is zero. 


The nonlinear equations of motion (6.17) are solved iteratively to 

determine displacements and accelerations at t+At beginning with a known 

solution for time t. Suppose that ’m' such iterations have been 

performed and the current estimate of the total nodal displacement field 

is ^ +A1 "{a } and the nodal acceleration field is ++A ^"{a }. 

m m 

If t+At {R ffl } computed using (6.17) vanishes for these displacements arid 
accelerations, the solution has converged. Computations for the next 
time step are begun. If {R^} does not vanish, a generalized force 
imbalance exists at the nodes and is given by the terms in }. 

These are usually termed the residual loads; they must be eliminated by 
suitable changes of nodal displacements and accelerations during the m+1 


i teration. 
necessary 


to 


Let {AR^} represent the change of the residual 
| A i* 

remove {R^}; then the following equation is 


I oads 
to be 


solved 
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++At m } + t+4t {AR } = {0} 
m m 

The increment of the residual load vector is simply 

summing the increment of each vector in (6.17). 

++4t {AR } = t+4t {AP ) - ++4t {AF > - [Ml t+4t {Aa } 
m m m m 

•in which {P} and {F} are implicit functions of 

displacements. Corrections in the nodal accelerations, {Aa }. are 

m 

coupled indirectly to the nodal displacement corrections through the 
integration operator. The finite size corrections implied in (6.19) are 
also nonlinear functions of the spatial variables. By assuming that an 
approximate correction can be obtained using the differential of {R}, 

++At {AR } » ++A+ {dR }, (6.20) 

m m 


determined by 

(6.19) 
the nodal 


a one term Taylor series about the conf iguration {a^} yields a set of 
linearized correction terms for use In (6.19) 


dP 

++A+ {AP } * ++A+ {dP } = [r-2] {a } = [K,l {a } 
m m L da J m [ m m 


ttAt 


{AF } - ++A+ {dF } = [ 7 --] (a } = [K T ]{a } 
m m L da J m T m 


dF 


( 6 . 21 ) 

( 6 . 22 ) 


Vector {a^} denotes an increment of nodal displacement. The symbolic 
differentiations inside the [ ] lead to the initial load stiffness, 
[K 1, and the conventional tangent stiffness, IK ]. These two 

matrices are the Jacobians of the functions {P} and {F}, which are 

nonlinear in the nodal displacements. Explicit forms of [K^l for the T, 
L. and U. L. formulations are provided in following sections. The r m f 
right subscript outside the [ ] denotes the nodal displacements at which 

00 

the stiffnesses are evaluated. The influence of {Aa } on {AR } is 

m m 
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simply a function of the integration operator and need not be considered 

in (6.21) or (6.22). As discussed in Section 5.4, the [K^ 1 matrix 

arises due to non-conservative loading. It is generally non-symmetr i c 

and thus iqnored in forming {AR }. Corrections in the effective 

m 

loading, {P}, due to geometry changes are incorporated by occasionally 
recomputing {P} using the current geometry. Symmetric forms of t K^_ ] may 
be derived for both T. L. and U. L. formulations as shown in subsequent 
sections. 

Substitution of (6.22) and (6.20) into (6.18) and then into (6.17) 
yields the incremental -iterative equations of motion 


[ K ] { a > = ++A+ {P } - ++A+ {F> " [Ml ++A+ {a } 
r m m m m 


(6.23) 


which can be solved for {a^} using linear equation solving techniques. 

The integration operator enters (6.23) through the {a^} term. Each 

*t*h 

vector on the right hand side of (6.23) is computed for the m 
iteration estimate of the nodal displacements at t+At using (6.8), 
(6.10), (6.14) at the element level and (6.16) at the structure level. 
The improved estimate of nodal displacements for t+At is 


t+At. , t+At r , . x 

m+l m m 


(6.24) 


** 

Correspond ing accelerations are computed using the 
integration operator. The right side of (6.23) is then evaluated for 
the improved displacements and accelerations (the R.H.S. is 
simply t+At {R}). If this vanishes, the solution has converged; 
otherwise a new correction (a m ) is computed and the process repeated. 

Simplifications in solving (6.23) involve the use of various forms 
of the Newton-Raphson procedure. It is not necessary to recompute and 
triangulate [K^. ] during each iteration. Any previously computed IK^l, 
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including the linear elastic [K 1 , may be employed. The convergence rate 
of the iterative procedure will be affected accordingly. Over- 

relaxation techniques used in static analysis may also be used here to 
enhance the convergence character i sties. 

In a multilevel substructure model, (6.23) is applied only at the 
highest level structure which contains the synthesized linear substruc- 
tures and nonlinear finite elements. The effects of substructur ing are 
reflected in each term of (6.23). Modal synthesis techniques discussed 
in previous chapters provide procedures to formulate [Ml and [K] at the 
outset of the transient response computation. The external dynamic 
loading is generally independent of the nodal displacements 

(conservative loading) but not independent of time. Thus, if time 

dependent loads are applied inside a synthesized substructure, the 

equivalent dynamic loads on the nodes (and modes) remaining after 
synthesis must be recomputed each time step. A simplification occurs 
when the spatial loading pattern applied inside the substructure is 
synthesized once, then a specified time function is applied to the 
equivalent loads to yield the variation with time. The internal 

A -j*' 

resisting forces, {F}, are the sum of contributions from the 

individual nonlinear elements and from the synthesized substructures 

that appear in the highest level structure. Nonlinear element 

contributions are given by (6.8); synthesized substructure contributions 

are obtained from the product of their stiffness matrix with 

the m"*"* 1 estimate of the total nodal displacements *" +A "*"{a }. This 

m 

computation requires the retrieval of each synthesized stiffness matrix 
from secondary storage during every equilibrium iteration. The major 
computational savings derive from the greatly reduced size 
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Initial Calculations 


1. For each linear substructure (repeat at each level for nested 
substructures) : 

a. Form C«] using standard techniques, 

b. Form CM] as consistent or lumped. 

c. Reduce Ck] and [X] using any of the dynamic reduction techniques 
described in Chapter 3. 

2. Form the synthesized stiffness, CK1, and mass, CM], fo r the highest 
level structure. 

3. Set initial d isp lacements °{a}, and velocities, °{a}. 

4. Compute initial accelerations, {a}, from equi I ibrium equation: 

Cm] °{a} + DO °{a} = °{P}. 

5. Define constants for the specific integration operator, e.g., 
Wilson-0, Newmark-3. Denote these constants a., a., a., etc. 

i j k . 

6. Compute the contribution of the mass to the effective stiffness for 

the highest level structure: CK] = a. CM]. 

7. Triangulate effective stiffness, CK p ] = CK] + CK]> of the highest 

t j 

level structure, using Choleski decomposition: CK^] = CGCL] . 


Table 6.1 — Procedure for Transient Analysis with Substructur i ng 
(Damping Negiected) 
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For Each Time Step t -> t+At : 

1 . If specified by the user, update the tangent stiffness, of the 

highest level structure. Only elements that are currently nonlinear 
are updated. Condensed substructure stiffnesses are re-used. 
Triangulate the new CK^.] = CLUCL]^. 

2. Compute the effective load increment vector. Place it in {R}. 

++a+ {r} = ++A+ {p}_ + { P } _ [X](a. + {a} + a . + {a} ). 

^ J 

This expression is derived by subtracting (6.17) evaluated at time t 
from the same equations at time t+At as follows: 

Cm: ++A+ {a} + ++a+ {f} - Cm] + {a} - + {f} = ++A+ {p} - + {p}. 

Now substitute CK-p]{Aa} = ^ +A "*"{F} - + {F} and ^ +A "*"{a} = 

”f° O *1” DO 

a.{Aa} + a. {a} + a, {a} from the integration operator to yield 
1 J K 

( CK t : + a. C m: ){Aa} = ++A+ {P} - + {P} + Cm:( a. + {a> + a. + {a> ) 

The right hand side is simply ^ +A "*"{R}. It is not necessary to include 
any remaining residual load, ”*"{R}, in the effective load vector for 
the new step if equilibrium iterations will be performed. Otherwise 
+ {R} should be added to the ++A "*"{R} derived above. 

3. Solve for the displacement increment vector CGCO (a} = {R}. 

4. Proceed to step 5 if no equilibrium iterations are specified to correct 
the acceleration or to eliminate residual loads due to nonlinear 
response in the spatial variables. 


Table 6.1 -- (continued) 
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4. (continued) Otherwise, begin equilibrium iterations: 


a. Set m = I 

b. Compute improved estimates of the tota I nodal displacements and 
accelerations given values at time t and the most recent estimates 
of change over the step. 


++A+ {a } = + {a } + (a }; 
m m 


++A+ {a } = a. ++A+ {a } - a. + {a} - a. + {a}. 
m i m j K 


c. 


For the improved total displacements and accelerations, compute the 

total residual load vectort +A tR } = ^ +A ~^{P } - "*" +A ^"{F } - CM] ^ A+ {a m }. 
m m m m 


d. Solve for a correction to the displacement change over the step: 

[LXa T {Aa} = ++A+ {R m }- 

e. Update the estimated change in nodal displacements for the time 


step: {a .} = {a> + {Aa}. 

m+l m 

f. Perform checks on convergence using {R m }, ^Aa}, (a m ). 

converged, jump to step 5. If not converged and iteration number 

m is less than the maximum allowed, increment m and go to g; other- 
wise, terminate analysis as a nonconvergent system exists. 

g. If specified by the user, update [K.j.1 and triangulate. Go to b. 

5. Compute new accelerations, velocities, and displacements for t+At to 
serve as initial conditions for the next step: 

++A+ {a} = a. + {a} + a.^ja} + a, + {a}. 

1 J K 

t+Atj-*-, t,*, , t/\ , t+At 

{a} = {a} + a. {a} + a {a}. 

I m 

t+At { i +r \ , r -i 
{a} = (a) + {a}. 


Stresses, strains, reactions, etc. can be computed for the highest level 
structure. Similar results may be recovered for condensed substructures. 


Table 6.1 


(conti nued) 



229 


of [K_] possible with substructuring. The frequent tr iangu I at ions 
of [K^l required to enhance convergence each consume much less time than 
does the reduction of a IK^] representing an equivalent unsubstructured 
model . 

Table 6.1 provides a detailed flow of the transient response compu- 
tation for a multilevel substructured model based on an implicit 
integration operator. 

6.4 Total Laqranqian Stiffness 

Specific forms of the elemental stiffnesses for the T. t. approach 
are described in this section. Once element stiffnesses are generated, 
a substructure stiffness is assembled using standard, well documented 
techniques [6.121. 

The approach adopted here derives from the' work of Nayak [6.21, 
McKay [6.131, and Dodds [6.141. As shown in (6.8), the element nodal 
force vector due to internal stress and strain is given by 

{F e } = / ,. *[B]*{cj}*dV (6.25) 

# v 0 

in which the above integral may be evaluated at any time t. With the 
configuration at time t known, an approximation for {F } at t+At is 
g i ven by 

++A V e } « f {F e } + d + {F 6 } = + {F®} + + [K®]d{a} (6.26) 

o 1 

in which d{a} represents the change In element nodal displacements over 
At. 

For the T. L. approach, the (*) configuration in (6.25) is taken to 
be that at time 0. The integrand is evaluated at time t but referred to 
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the configuration at time 0 using the results derived in Chapter 5, 
specifically (5.81). After substituting into (6.25) the differential, 
internal force vector may be written as 

d + {F e } = d[ / v *[B] T *{S}°dV] (6.27) 

o 

Using the product rule, differentiation under the integral sign yields 

d + {F S ’} = / .. + [B ] T d + {S> t d + IB ] T + {S}°dV (6.28) 

J v o o o o 

o 

From (5.100), the increment of 2nd P-K stress is given by 

d + {S} = + ED T Id t {s} (6.29) 

o o T o 

with the tangent modulus matrix provided by (5.101). The differential 
of the Green strain is given by (5.55) as 

d + {e} = + [B]d{a} (6.30) 

o o 

t 

in which the form for q [B1 is repeated in Table 6.2. 

The differential of the ^fB] matrix in (6.28) with respect to the 

o 

element nodal displacements has been derived in a convenient form by 
Nayak [6.21 and Dodds 16.141. The *[BI matrix is split Into two 
components, the first being independent of the nodal displacements which 
vanishes under differentiation with respect to the nodal 

displacements. Differentiation of the nonlinear terms yields an 
expression of the form 

d + [BJ t {S} = + [G1 T+ (M1 + [G1 d{a> (6.31) 

o o o o o 

following rather lengthy but stra ightforward differentiation of the 
summations. Complete details of the process are given by Dodds 
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1 ” 

[6.14]. Matrix [G1 simply contains derivatives of the element shape 
functions with respect to the initial, °{x}, coordinate system. 

4 * 

The [Ml matrix contains a regular arrangement of the 2nd P-K stress 
components, ^{S}. Details of the [G ] and [ M ] matrices for the 2-D case 
are given in Table 6.2. 

The complete element tangent stiffness is obtained by combining the 
effects of (6.28-6.31) as 

f {AF e } = j’tK-.KAa} = [X, ] t + [K ]]{Aa} (6.32) 

o T L o NL o a J 

in which finite increments of the nodal displacements, {Aa}, are used 
during computation to replace the differential values. 

The initial stress stiffness, [K ] , is given by (6.31). The 
nonlinear stiffness, [K^], is given by 

+ [K... 1 = / „ + IBI T tED T ]'tfBl°dV e (6.33) 

o NL i V o o T o 
o 

and contains the usual linear stiffness, in addition to displacement 
dependent contributions, i.e., if the nodal displacements are 

zero, q [B 1 simplifies to the standard linear strain-displacement rela- 
tion. 

6.5 Updated Laqranqian Stiffness 

The procedures adopted in the previous section also enable the 
derivation of stiffness matrices for the U. L. formulation. The 

resulting stiffness matrices for U. L. have a form very similar to that 
for T. L. which simplifies the associated programming when both 
approaches are implemented. 


232 


The Green strain increments in 2-D are: 


^{6e} = ^]B]{6a}, 


/{£} = *[B] d{a} 


{e} T =1 + e ; + e ; + Y ~| {6a} T = I 5u ; 6v I ; variation of total 

Lox'oy’ ° Y xy I L I disp | acemen+s 


t t , . -r t 2 , t 2 -| 

e = u + 0.5 l u + v ]; 
ox ox ox ox 

t t , n , r t 2 + + v 2 "]; y = °y 

£=v+0.5Lu o y ' 7 7 

o y o y o y 7 

t t t t t t t 

Y = U + V + U U + V V 

oxy oy ox oxoy oxoy 


t 


~t 
3 u 


x = x; u - „ 

ox „o 
3 x 


A virtual or differential change of the strain £ is given by 

U /\ 

4- -4* "J* “4- ~j~ ■j* ( 

6c =5u + u6u + v6v with other terms derived similarly, 

ox ox oxox oxox 


IB] = 


(I + + u )i- 
° x 3x 


t 3 ' 
u — 

° Y 3y 


t 3 


« t 3 
I v — 
ox,.. 

I 3x 


I (| + + v ) 1. 

: oy 3y 


(I + + u ) — 
ox 3y 


(I - + v ) i- 
° y 3x 


+ v i- 
° x 3y 


Table 6.2 — Summary of Total Lagrangian Matrices 
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Substituting shape functions into the above, the corresponding matrices for 
the Jth node are: 


- 


t 3N J 
(I + u )— 

° x 3x 


t 

o u y 


3y 


t 3N J 
u — — 

0 y 8x 


t 3n j 

( i + u ) — 
ox ~ 
3y 


t 3N J 

v 

o X ~ 

3x 


t 3N J 

( I + V ) — 

° y 3y 


t 3N J 
(I + v ) — 

° y 3x 

° V>< 3y 


= C V ; V ; 


. + r n 1 
; o G 11 


y j ^ T 


CM] = 



; 3n j 

3x 

1 3y 

i 

1 

0 

: * 
i 

s 

i s 

xx 

i xy 

S 

i 

1 s 

xy 

; yy 

0 

i 

0 

i 

! 0 

. 

i 


-4. 


1 3N . 

i ± 

1 3x 
i 

0 I 
_ i 

0 


9N , 


3y 


0 

0 


± _ 
I 

I 

I 

■ r 


s y s 

xx I xy 


S IS 
1 xy , yy 

i 1 


i J 


+ S„. 
o 'J 


Table 6.2 — (continued) 
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The basic equation of the U. L. formulation defines the internal 
virtual work of an element at time t+At relative to the configuration at 
time t as 


t+At r ,.,e 

int 


SWT ^ = -/ v t+ ^ + {5e}' T T'{S}'dV 


i t+At r^it. 

t 


(6.34) 


Defining a suitable differential operator matrix, ++A ^[B1, to yield vir- 
tual strains from virtual nodal displacements, (6.34) may be written in 
the form of (6.8) as 

t+A+ SW® n+ = ~{<Sa e } T / W ++A t [B 3 1 ++A t(S} + dV = -{5a e } T ++A+ {F e } (6.35) 


.V f 


in which ^ +A ^{p e } In (6.35) has the identical meaning as the same 
quantity in (6.26). 

The element resisting force vector at t+At is expressed in terms of 
a known value for the configuration at time t plus a differential change 
over At as 


t+At rr e. t r e. , ,t r _e, 
{F- } ** {F } + d {F } 


(6.36) 


In which 


dJ{F®} = JtK T ]d{a} (6.37) 

The difference between (6.36-37) and (6.26) is simply the reference 

*t" © 

configuration of the tangent stiffness. The incremental force, d {F }, 
and nodal displacement Increment, d{a}, have the identical physical 
Interpretations in both T. L. and U. L. formulations. 

•j" 

The specific form of is obtained by evaluating the 

differential of (6.35) with respect to the displacements at the value of 
nodal displacements correspond ing to time t. Symbolically, 
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d f { F e > = / + ^ [B l T d^{S> + d*[B] T *{S} + dV (6.38) 

From Chapter 5, the following relationships are used to expand (6.38) 

^{S} = "^{a} (the Cauchy stress) (6.39) 

d^{S} = + {S} = ^(D T ld^{e} (6.40) 

The tangent modulus matrix, referred to the current (time t) configura- 

-f* 

tion, [ D-p 1 , was derived in Chapter 5. The differential strain, 

measured from the current configuration, is 

d*{e} « ^IBId{a} (6.41) 

in which ^[Bl is simply for a zero increment of displacement 

(see Table 6.3). Notice that the resulting form is the conventional 
linear £B1 matrix containing derivatives taken with respect to total 
displacements at time t. Using (6.40) and (6.41), the first integrand 

•j* 

term of (6.38) becomes the "nonl inear" stiffness matrix ,[K.„ ] and is 

t NL 

given by 

t lK NL ] = ft t tB L lT JlD T 1 + EB L 1+dV (6.42) 

in which the "L” subscript is used to denote the simple, linear form 
taken by [ B ] . 

To obtain the second integrand term of (6.38), an expression 
for d^JBI is required. This is found by differentiating with 

respect to the nodal displacements, then evaluating the terms for a zero 
displacement increment. A comparison of ^[Bl in Table 6.3 

-f 

with I B I in Table 6.2 shows a very similar form. The same operations 
that yield d^tBJ, readily yield d^[B] as 
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The Green strain increments in 2-D are: 

t + At We) = t + At [B]Ka}> d t+A+ {e} , t + At [B] d{a) 

++ >> T ’ t {e>T ’ |_t e x ' Ty : t^xy"~ | 

{ <Sa ]- ^ = 6 u ; <5,v I ; variation of incremental displacements that 

— ' ' I produce configuration at t+At. 


. u- 

t X 

1 1 

in 

o 

+ 

CM 1 X 
+- 

+ ,V- ']; 
t X ’ 

x = °x + u ; + u- - 

9^_u 

9x 

t V y 

+ 0.5C 

2 

, u- 

+ y 

+ + V- ']; 
t y 

> 

+ 

>- 

0 

ii 

1 



t Y xy - t u y + t v x + t u x t u y + t v x t v y 

A virtual or differential change of the incremental strain is given by 

6 e = j u- + , u~ 6 L u- + , v- S.v- ; with other terms derived similarly, 
t x tx txtx txtx 
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Substituting shape functions into the forgoing, the corresponding matrices 
for the Jth node are: 


++A |Cb j ] = 


( I + ,u-) 
t x 


3N 


3x 


+ u y 


, u- 
t y 


3y 

3x 


(' + t u 5 ) 


3N 


3y 


t V x 


3N 


3x 


3N 

( , + v ) 

Y 3y 

3N 


, v- 
t x 


3N 


3y 


CB] 


3N 


Jl 

3x 


3N 


3y 


3N 


3y 


3N 


3x 


It 


;ct] T 


Cm] = 


CO 

+- 4- 

1 

> 

+ r 2 . 
t G ’ 

• 

® © • 

t p n 
; t G 

1 



i 3Nj 

i 

1 0 

i 

1 

0 


3x 

1 

. _L _ 

3y 

1 

1 

_ L _ 

, 


0 

1 

1 

0 

1 3N . 
i J_ 

I 3n j 


_ 

1 

1 


\ 3x 

i 3y 

0 xx I 

T 

xy 

i 

L 

0 1 
1 

0 


T xy 1 

a 

yy 

i 

i 

0 1 
1 

0 


0 1 
I 

0 

i 

i 

a 1 

XX | 

T xy 


c 

i 

0 

1 

J 

T 1 

xy , 

Q 

■< 

< 



a. . = ,a. . 
|J + i J 


Table 5.3 — (continued) 
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d*[B] |{S} = f f lG] J j[M] ^IG] d{a) (6.43) 

in which the details of 1[G1 and "fjM] are shown in Table 6.3. 

t t 

Matrix ^[M] contains the Cauchy stresses at time t arranged in the iden- 
tical form that the 2nd P-K stresses appear in ^[Ml. Similarly, the 
t t 

form of ^[Gl is identical to that of q [G 1 except that derivatives are 
taken with respect to the displacements at time t rather than time 0. 
The integral of the product shown in (6.43) is again termed an initial 
stress stiffness. 

The complete element tangent stiffness is obtained by combining 
(6.42) and (6.43) to yield 

f {AF e } = *[K T l{Aa} = [ + CK Nl 1 + ] {Aa> (6.44) 

6.6 Comparison of Formulations 

Details of a transient solution procedure based upon an implicit 
integration scheme have been described in this chapter. An implicit 
scheme is recommended over an explicit scheme for three reasons. First, 
the procedures for static linear and nonlinear analysis are recovered 
from the implicit scheme by simply omitting the mass matrix; the 
degenerate case of an explicit scheme does not yield a formulation 
suitable for static analysis. The capability to perform both static and 
transient analyses with the same software is particularly attractive to 
engineers since a static analysis invariably precedes a dynamic 
analysis. The second reason to select an implicit scheme is the more 

general class of structures that may be analyzed. An Implicit scheme 
may be used to compute the details of localized wave propagation under 
very high velocity Impact as well as the vibration characteristics of a 
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massive structure subjected to time dependent loads. The explicit 
scheme may be computationally more efficient for wave propagation 
studies but not necessarily more efficient than an implicit approach for 
the latter class of problems. Thirdly, the equilibrium imbalance due to 
nonlinear response in the spatial variables can be "iterated out" to 
within a specified tolerance using the implicit scheme. Dynamic 
equilibrium can thus be assured at the end of each time increment. 

Complete details of the element tangent stiffness matrices have been 

presented in the 2-D case for both the T. L. and U. L. formulations. 

For both approaches, stiffness matrices have the identical symbolic 

form, which is conveniently expressed as the sum of an initial stress 

stiffness, [K 1, and a nonlinear stiffness, [K.„ ]. Given the common 
a NL 

rate of work per unit mass expression from which each formulation is 
derived, it is expected that identical results for each solution would 
be obtained (provided the full stress rate transformations described in 
Chapter 5 are utilized). Some computational evidence [6.10, 6.111 does 
demonstrate the agreement of overall load-deflection curves for static 
and dynamic solutions which include plasticity effects. Unfortunately, 
no comparisons are provided of the internal stress distributions, which 
must be identical if equivalence of the formulations is to be 
demonstrated. Global agreement of load-deflection does not imply 
identical internal stress distributions. 

The computational efficiency of each approach is addressed at two 
levels; namely, the effort required to compute the element matrices, and 
the number of iterations required to attain the equilibrium 
configuration at time t+At, assuming At is the same in both 


formu I at ions 



240 


!n each formulation the construction of IK ] requires the same 
number of operations if the expense of coordinate updating in U. L. is 
neglected. In contrast, construction of the nonlinear stiffness, 
[K nl ], in U. L. requires the same number of "operations as required for 
a conventional linear stiffness. The T. L. nonlinear stiffness requires 
.more operations due to the absence of zeroes in the IB] matrix (refer to 
Table 6.2). The increase in the operation count depends on the element 
type but easily exceeds by a factor of three the number of operations 
required to form IK^] in U. L. Evaluation of the internal resisting 

force vector during iterations requires the identical number of 

operations in each case; the U. L. matrix has the identical form 

of ^[Bl in T. L. (no zero terms). The reduced operation count of U. L. 
compared to T. L. during element stiffness generation is an important 
consideration only when a 1 1 matr i ces associated with the analysis reside 
in memory during execution. Once the swapping of element stiffnesses to 
and from disk begins, it would appear that the I/O overhead overwhelms 
any advantage of one formulation over the other. Furthermore, the 

timing results shown in Chapter 2 for several large linear and nonlinear 

analyses reveal that element stiffness generation times do not represent 
a major portion of the total solution time. Thus, a reduction of 
stiffness generation time obtained by the selection of U.L over T. L. 
does not yield an equal percentage reduction In total analysis time. 

Perhaps the most Important efficiency comparison between the two 
formulations is the number of equilibrium iterations required for 
convergence during each time step. Using a common program that has both 
U. L. and T. L. capabilities, a structure could be analyzed with each 
formulation for identical time steps and with an identical convergence 
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criterion based on total force equilibrium. The results of such a 
comparison would demonstrate any inherent computational advantage of one 
formulation over the other, a least for a particular problem. 
Unfortunately, no such comparisons have been found in the open 
I tterature. A systematic study of several problem classes using this 
.approach is necessary to form the basis for any general efficiency 
statements. 

The distorting of elements has led to numerical problems in some 
U. L. finite element analyses. For example, an element that is square 
at t=0 may become a badly distorted quadri lateral at time t+At. Element 
response character i st ics are known to degenerate rapidly as aspect 
ratios increase; this will influence the response in an U. L. approach 
as shape function derivatives are dependent on the deformed shape at 
time t. In T. L«, all such effects are incorporated in the pre- 

f 

multiplier terms of shape function derivatives in the ^[Bl matrix. The 
shape function derivatives in T. L. are always computed relative to the 
configuration at time t=0 and thus remain well-behaved if the aspect 
ratio at t=0 is acceptable. The interesting case that demonstrates the 
problem with U. L. is an 8-node, 2-0 i soparametr i c element which Is 
square at t=0 but whtch at time t has the mid-side node displaced toward 
a correspond ing corner node. Any movement of the mid-side node toward 
the corner produces a singular point In the correspond ing shape function 
when the derivative is evaluated with respect to the current element 
shape. No such singularities are Introduced in shape function 
derivatives in T. L. if the mid-side node is properly positioned at t=0. 

A final comparison of the two formulations considers their applica- 
bility in a substructured modeling and solution procedure. At all times 
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t, the stiffness matrix of a linear substructure and its synthesized 
form remain unchanged from that at time t=0. The configuration of the 
structure at t=0 provides the basis to which all strains and stresses 
are referred. The linear, strain-displacement relations and constitu- 
tive matrix that is independent of nodal displacements at time t are 
utilized. Consider the common nodes along the boundary of a nonltnear 
region that is shared by an adjacent I inear substructure. In an U. L. 
formulation, the coordinates of these common nodes are updated at time t 
to reflect the incremental displacements over the previous increment 
At. The tangent stiffness of the nonlinear region is regenerated to 
reflect the new geometry. Thus, there exists a discontinuity of nodal 
positions along the boundary. If the substructure stiffness is 
recomputed and synthesized for the updated geometry along the boundary 
(and the consequent Internal repositioning of nodes) all the advantages 
of a substructured model are lost. The alternative is to assume that 
nodes along the interface are sufficiently remote from any effects 
causing nonlinear behavior that the response In nonlinear elements on 
the boundary Is actually linear, i.e., displacements and strains are 
infinitesimal and the matrices involved at the element level revert to 
the linear analysis form. If such an assumption is actually verified In 
the analysis, the integrity of the results is assured. 

In a substructured T. L. solution, discontinuities of nodal 

positions along a boundary do not occur since nodal coordinates are 

never updated for the displacement increments. However, the effects of 

deformation and rotation in the nonlinear elements adjacent to the 

boundary are Implicitly incorporated in the tangent stiffness through 

1 * 

the displacement dependent terms in [BI. Thus, the same difficulties 
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described above for U. L. appear to arise for T. L. However, the 
maintenance of external geometric continuity is more appealing than 
allowing the discontinuities to develop. It could possibly be argued 
that the T. L. formulation would introduce less disturbance in the 
strain field across the boundary than II. L., but there exists no 
computational evidence to collaborate this speculation. 

In view of the considerations described above, a Total Lagrangian 
approach is recommended to form the basis of a general software system 
having broad applicability. If computational evidence in support of 
LI. L. as a more efficient approach should be forthcoming, only minor 
changes in software to support a To L. formulation are required. In 
fact, nearly all U. L. dependent computations can be isolated within 
element dependent routines. The reverse situation of implementing a 

T. Lc approach in a system designed only for U. L. Is not nearly as 
simple. 
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CHAPTER 7 


USER INTERFACE — INPUT DESIGN 


7.1 General 

The most popular approach to user communication with structural 
analysis software is the problem oriented language (POL). Virtually all 
successful software systems use the POL approach, either by Initial 
design or by the use of pre- processors to translate POL input Into fixed 
format, card Images. The POL approach provides the user with greater 
flexibility by placing him in control of the process rather than forcing 
him to conform to rigid formats and sequences. The self-documenting na- 
ture of the input reduces the need for reference to manuals and provides 
a concise description of the structural model for other analysts. The 
POL Is essential for interactive processing where error recovery is 
often necessary. 

Dynamic analysis with multilevel substructuring will be Implemented 
as an extension of the present POLO-FINITE structural mechanics system. 
The philosophy established during the development of POLO-FINITE was to 
maintain as much Independence as possible among the various components 
of a complete structural model. These components include nonlinear 
material model specification, geometric definition of structures, 
parameters controlling the nonlinear algorithms, and requests for com- 
putation and output. The primary reasons for choosing this approach are 
to provide maximum flexibility in using condensed substructures as ele- 
ments In higher level structures and to minimize the effect of changes 
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In the structural model throughout the analysis-design sequence. 

Wherever possible, this philosophy is maintained in the extension 
to dynamics. However, one area exists in which dynamic solution 
parameters must be tied directly to the geometric definition of a sub- 
structure. This is the frequency analysis of a substructure that Is to 
be condensed by modal synthesis. Since economical frequency analysis 
depends upon the type of structure, the number of eigenpairs required, 
and the solution method. It Is not appropriate to select just one solu- 
tion algorithm for all substructures in a complex model. Various sub- 
structures will have differing character i sties and may require an une- 
qual number of retained modes for condensation. It Is also possible 
that one substructure could be condensed two or more times In differing 
ways, with varying geometric and generalized DOF, for use In separate, 
higher level structures. Thus, It Is necessary to tie the selection of 
the etgenproblem solution method to the structure definition. 

The capabilities selected for general purpose dynamic analysis, 
along with the various options and parameters that control the solution, 
must be defined accurately and unambiguously by the POL. Section 7.2 
presents an explanation of the capabilities to be Incorporated Into 
POLO-FINITE. Section 7.3 lists the command structure for dynamics. Ex- 
amples Illustrating the use of these commands are given in section 7.4. 
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7.2 Description ±M EQL 
Structure and E I ement Mass 

The mass of a structure can be divided Into two parts; primary and 
secondary. Primary mass Is the mass of the load-carrying components 
(elements) of the structure. Its definition Is easily added to the 
specification of an element through two new element properties- The 
first defines the type of mass formulation; LUMPED or CONSISTENT. The 
second is the DENSITY of the material of which the element Is composed. 
The element mass matrix can then be formed using existing shape func- 
tions. Assembly of primary mass for a structure will follow a procedure 
identical to that used In stiffness assembly. The current FINITE system 
accepts up to thirty DOF at each node of an element. These are the 
displacement DOF; U, V, and W, plus their first and second derivatives; 
UX, VX, WX, UY» etc. Depending upon the particular element formulation* 
it Is possible for mass to be assigned to any or al I of these DOF. 

Secondary mass Is the mass of non- load-carrying components, such as 
concentrated and distributed live-loads, that are supported by the 
structure. Secondary mass Is defined In a manner similar to the ap- 
plication of gravity loads. The secondary mass Is resolved Into equiva- 
lent nodal masses via the appropriate shape functions. The result will 
always be a lumped mass matrix which Is added to the primary mass of the 


structure. 
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There are three types of secondary mass: nodal, element, and pat- 
tern. Nodal mass Is mass that Is concentrated at a structure node. 
Element mass Is concentrated or distributed on the surface of and ele- 
ment. As with primary mass, secondary mass may be assigned to any of 
the thirty nodal DOF. The pattern mass Is provided as a convenience to 
the user. It enables the definition of secondary mass In terms of a 
previously defined loading condition, usually gravity loading. The user 
need specify only the name of the loading condition to be used as the 
pattern and a value for the acceleration of gravity to support the ap- 
propriate conversion from force to mass. 

Substructure loads can be defined at the lower levels and then ap- 
plied selectively at the higher levels of the structure hierarchy 
through the "EXTERNAL ELEMENT LOADS” facility. It would be advantageous 
to have a similar capability with respect to substructure mass. The 
analyst may wish to use several copies of a particular substructure, 
each with a different mass distribution (described by secondary mass In- 
put). This analogy to substructure loading Implies the need for an "EX- 
TERNAL ELEMENT MASS" facility. This Is not possible in dynamic analysis 
since the change In the mass of a substructure changes the natural fre- 
quencies and mode shapes. Physically distinct substructures exist when 
the mass distribution varies. Each of these distinct substructures must 
be uniquely defined at the lower levels of the hierarchy- The mass is 
then automatically carried through the hierarchy via the condensation or 
synthesis process. 
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The commands for computation (assembly) and output of the mass 
matrix for a structure or stand-alone element fellow directly from those 
for the stiffness matrix. 

JS ±o i£±ii £a .Damping 

Since the dynamic reduction process recommended here does not in- 
clude the substructure damping matrix, damping is defined only for the 
highest level structure. Two methods are available for defining struc- 
tural damping: modal and Rayleigh. Definition of modal damping re- 
quires input of the modal damping ratio for each vibration mode under 
consideration. Modal damping Is applicable only to transient analysis 
by mode superposition. Rayleigh damping involves the definition of two 
damping ratios at two selected frequencies! the frequencies need not be 
eigenvalues of the structure. Rayleigh damping Is applicable to tran- 
sient analysis by either mode superposition or time-history integration. 
Use of Rayleigh damping requires that a frequency analysis be performed 
in order to compute the modal damping ratios for mode superposition or 
to explicitly form the damping matrix for time-history integration. 

Depending upon the method used to define damping, either the 
damping matrix or modal ratios can be output for the structure. 

Frequency AnaJyS-Ls 

As previously mentioned, the parameters controlling the frequency 
analysis (computation of natural frequencies and mode shapes) must be 
defined explicitly for each structure for which the analysis is to be 
performed. No default analysis method is adopted. The syntax for 
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specification of the solution method Is similar to that for a nonlinear 
material. That Is, the TYPE of solution procedure Is Identified fol- 
lowed by a listing of PROPERTIES unique to that type. When appropriate 
the range of frequencies and the maximum number of modes to consider are 
specified at this time rather than In the computation request. Solution 
method properties can be changed via analysis restart and the ACCESS 
STRUCTURE... sequence. If a substructure Is to be condensed by Guyan 
reduction, no frequency analysis specification Is required. 

In the request for computation the analyst may select a nonlinear 
dynamic loading and time step at which the frequency analysis Is to be 
performed. This allows the user to Interrupt a transient analysis after 
some nonlinear behaviour has occured and compute natural frequencies and 
mode shapes of the structure. Standard output Includes natural frequen- 
cies and modes shapes. 

Prior to a transient analysis by mode superposition, the user may 
examine the modal content of a particular dynamic loading condition- A 
special output request facilitates selection of the modes that par- 
ticipate In the dynamic response. After a frequency analysis the user 
can request output of MODAL LOADS for the loading condition. The fre- 
quency content of the loading can then be examined and the appropriate 
modes selected for superposition. 
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Dynamic R educti o n 

The procedure to request dynamic reduction paral lets that for 
static condensation. The reduction method is defined at the inter- 
mediate substructure level j i.e., the substructure with only one element 
of type: CONDENSED. Guyan reduction is the default method. Automatic 
selection of master nodes, in addition to the Interface nodes, is In- 
cluded by specifying the number of additional nodes to be retained- The 
f Ixed- i nterf ace method is invoked by specifying which substructure nor- 
mal modes to retain. The modes specified must be within the range com- 
puted in the frequency analysis of the lower- level substructure. The 
retained modes need not be consecutively numbered. As an alternative to 
using normal modes, user-supplied modes can be included In the synthesis 
process. These modes could be derived from an experimental analysis or 
some other source, such as low-order polynomials. Input data describing 
these modes must be included with the definition of the structure to be 
condensed. 

Dynamic reduction can be explicitly Invoked with a COMPUTE STIFF- 
NESS... or COMPUTE MASS... command for the intermediate level substruc- 
ture. Reduction Is performed automatically when required to satisfy a 
request at a higher level. 
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Initial Conditions 

Initial conditions can be defined for a ‘structure prior to tran- 
sient analysis. They define a starting solution, in terms of displace- 
ments and velocities, for the unconstrained physical DOF at time t = 0. 
For all other times the displacements and velocities from the previous 
time step are used automatically. 

The user can specify Initial conditions In one of two ways. First, 
he can define numerical values for each DOF with non-zero displacement 
or velocity. The default Initial conditions are zero displacement and 
velocity for all unconstrained DOF. The second method uses the static 
equilibrium configuration from a previous linear or nonlinear analysis. 
This method allows the structure to be released from some deflected ini- 
tial shape with zero Initial velocity. A dynamic loading can then be 
optional ly appl led. 

fly n mt-C. leading 

The dynamic loading function, P(x,y,z,t), Is defined such that It 
has a spatially-varying component, F(x,y,z), and a time-varying compo- 
nent, G(+): 


P(x,y,z) = F(x,y,z) * G(t), 


( 7 . 1 ) 


Simply stated, the pattern of the load Is fixed and Its magnitude 
changes with time. 
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The load pattern, F(x,y,z), can be described as either actual 
forces applied to the structure or as support accelerations. The former 
can be best defined as a static linear loading condition, while the 
latter requires an additional loading type: NODAL ACCELERATIONS. No 
special provisions are necessary for input of out-of- phase support ac- 
celerations. They can be recognized and handled automatically. 

The time-varying component, G(t), is combined with other loading 
data to form a dynamic loading condition. The G(t) vs. t relation may 
be harmonic. Impulsive, or general. The dynamic loading condition must 
also include the loading pattern. F(x.y.z), which is to be used. More 
than one static linear loading condition can be combined to form the 
complete pattern of the dynamic load. Other necessary Input Includes 
the values of time t at which displacements are to be computed (thus, 
defining the step size) and values of time t at which computed results 
are to be retained in the data base. This last item is Important 
because a transient analysis of any significant duration could result in 
more data than could be effectively stored. Also, It is likely that 
stresses and strains would be required at only a few of the many time 
points for which displacements are computed. 

:[nm3J3Ji± A n a lys is . 

Transient analysis yields the response of the structure, In terms 
of displacements and possibly velocities, when It is subjected to time- 
varying loading or support accelerations. Two approaches are available 
for performing transient analysis: mode superposition and time-history 
Integration. Mode superposition requires that a frequency analysis be 
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performed so the equations of motion can be uncoupled. This implies 
that an appropriate frequency analysis method must be selected prior to 
requesting the transient analysis. The resulting set of independent 
equations is easily solved using one of the Lagrange interpolation for- 
mulas. Time-history Integration is performed by any one of a number of 
explicit or Implicit operators. Specification of the transient analysis 
method is similar to that for frequency analysis; the TYPE of method is 
defined followed by the PROPERTIES list. 

The request for computation includes the loading condition, time 
steps, and optionally Initial conditions and a mode list. The mode list 
is used with mode superposition to specify which modes participate in 
the response. Results available for output include displacements, 
velocities, strains, and stresses. 




Anal.ysJ ? 


The analysis of shock spectrum response is currently restricted to 
linear structures. The shock spectrum Is input by defining the func- 
tional relationship between a spatial coordinate and a time coordinate. 
The spatial coordinate can be chosen as displacement, velocity, or ac- 
celeration, while the time coordinate can be either period or frequency. 
The user inputs discrete points from the spectrum and the remainder of 
the curve is constructed by linear Interpolation in four-way logarithmic 
coordinates. The direction of application of the shock is defined using 
direction cosines for the translational DOF (U, V, and W for 3-D struc- 
tures). The nodes at which the shock Is applied are also defined. 
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Prior to computing the spectral response, a frequency analysis of 
the structure must be performed. Spectral response quantities are com- 
puted only after the corresponding output request has been made. 
Results available for output Include spectra! displacements, spectral 
velocities, spectral strains, and spectral stresses. These quantities 
can be output on a mode by mode basis or In some combined form. Methods 
used to combine the modal quantities include SRSS (square root of the 
sum of the squares) and PEAKJSRSS (peak response mode plus SRSS of the 
remaining modes). PEAK__SRSS Is also known as the Nava! sum. As a 
measure of the portion of the total mass responding to the shock In each 
mode, the modal PART I C I PAT 1 0N_FACT0RS can also be output. 

Utl I ity Caima n ds 

The dynamic solution process can be monitored by Invoking the 
TRACE... command. Messages listing the currently executing module and 
elapsed CPU time are output at various checkpoints. 

To eliminate unnecessary data from the data bases, the DESTROY... 
command Is expanded to include results from frequency, transient, and 
spectral analyses. 
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7.3 POL Structure 

7.3.1 Symbolic Conventions of the Syntax 

The following Is a description of the conventions used In this sec- 
tion to Illustrate the FINITE command syntax. 

A descriptor is used to Identify the position and class of a data 
Item in a particular FINITE command line. The descriptor Is delimited 
by the characters "< >.” The command 

NUMBER OF NODES <!nteger> 

Implies that the word NODES Is to be followed by an Integer. An ap- 
propriate example Is: 

NUMBER OF NODES 100 

The following are definitions of the descriptors used within the 

POL. 


<!nteger> — a series of digits optionally preceded by a plus or 
minus sign. Examples are 121, +300, -8. 

<real> -- a representation of a floating point number In either 
decimal or exponential form. Real numbers must contain a 
decimal points and may be optionally signed. Examples are 
1.0, -3.5, 5.24E-08. 

<number> — either an Integer or a real number can be Input. The 
data Item Is converted to a real number. 
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<In+eger list> — a sequence of integers. The sequence may be 
listed explicitly or defined over a range of integers with a 
constant increment. The default increment is 1. Examples of 
integer lists are: 1, 2* 4, 5, 8, 11; 1-10; 2 TO 20 BY 2. 

<real llst> — a sequence of real numbers. Real lists have the 

same form as integer lists except that there is no default 
Increment. Examples are: 1.0, 1.5, 2.0, 3.0; 0.0-2. 5 BY 

0.25. 

<number list> -- either an integer i 1st or a real list is input. 
The data i s converted to rea I . 

< I abe I > — a series of letters and digits beginning with a letter. 
Examples ares PLANEFRAME, DEADLOAD10. 

<string> — any text enclosed within single or double quotes. An 
example is: '‘THIS IS A STRING". 

in some instances a description of the physical meaning of the data Item 

is added to the class in the syntax of a descriptor- This Is helpful in 

clarifying the use of the data Stem. For example a command of the form 

STRUCTURE structure names I abe! > 

Implies that the data item following the word STRUCTURE Is a label 
defining the name of the structure. 
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It Is not always necessary to completely spell out every word on a 
command line in order to have the command correctly translated. Many 
words can be abbreviated and these are Identified in the command syntax 

by underlining. The underlined portions of words Identify the minimum 
Input necessary for proper command translation. Descriptors are not un- 
derlined but are replaced by an item of the specified class when ap- 
plicable. If the command syntax has the form: 

NUMBER OF NODE S < Integer > 

the following is acceptable as input: 

NUM OF NODE 10 

When only one word from a group of words may be selected as input, 
the choices are listed one above the other and enclosed In braces, 
"{ }". The command syntax 

COMPUTE j STIFF NESS j 

1 displ acements! 

implies that any of the following commands are acceptable: 

COMPUTE STIFFNES 
COMPUTE DISPLACEMENTS 


COMPUTE DISPL 
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When an entire word or phrase in the command Is optional, it Is 
closed within parentheses. The command with the syntax 

NUMBER (OF) NODE S <integer> 

can be written as 


NUM NODES 100 


When more than one word from a group of words may be selected, 
group Is enclosed in brackets, 3 H . The command 


OUTPUT 


DISPL ACEMENTS 


STRAIN S 

STRESSES 


Implies that the user may request 
OUTPUT DISPL STRAINS 

Brackets and braces are combined to allow more flexibility 
designing commands. The command syntax 


<1 nteger> 



Impl fes that the user may enter data of the forms 


1 X OoO Y 0.0 Z 5.0 


en- 


the 


2 X 1.0 


Z 5.0 
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Continuation of an input line onto a second physical line is accom- 
plished by placing a comma at the end of the line to be continued. 

Comments may be placed in the data by placing a C in column 1 and a 
blank In column 2 of the comment line. 

One method for line termination is to place a dollar sign "$" on 
the line. Space on the line following the H $" Is Ignored by the trans- 
lator and and may be used for comments. 


7.3.2 Command .Syntax 


fit Maa s , 


Example of the command to specify primary mass: 

ELEMENT 1 TYPE CSTR 1 ANGLE CONSISTENT E 1 NU .3 DENSITY .00074 

Example of the commands to specify secondary mass: 

MASS 

NODAL 

2 U V W 20.0 THETAX THETAY 5.0 
ELEMENT MASS FOR TYPE PLANEFRAME 

3 LINEAR U V W FRACTIONAL LA 0.25 LB 0.75 WA 3.0 WB 8.0 

1 CONCENTRATED UVW L3.6M5.0 

2 CONCENTRATED THETAZ L 3.6 M 3.0 
USE LOADING DEADJLOAD G 386.4 

Assembly command: 

COMPUTE MASS (FOR) ( STRUCT URE ) < I abe I > 

\ ELEM ENT I 


Output command: 


OUTPUT MASS (FOR) I STRUCT URE i < labe I > 

\ ELEM ENT j 
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&£ Damping 


Modal damping: 


DAMP ING MOD AL / RATIO S l [<mode I is+> <number>] 
\ PERCENT S f 


Rayleigh damping: 

DAMP ING RAYLEIGH [/ FREQU ENCIES^ <number> <number> 

( PERIOD ) 

( RATIO S ) <number> <number> 

I PERCENT S ) 

Output command: 

OUTPUT DAMP ING ( MATRIX ) ((FOR) STRUCT URE <labe!>) (,) 

{ RATIO S V 
t PERCENT S ) 

((FOR) MODES < Integer list» 

Units of seconds for time and rad/sec for frequencies will be required. 


Specification q± Frequency Analysis, 

Definition of the frequency analysis methods 

FREQ UENCY ANAL YS I S ( TYPE ) ( HQRI \ 

1 JACOBI / 

1 SUBSPACE I 

< NEWTON > 


PRO 
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Standard output request: 


OUTPUT (NONLINEAR) 


(NATURAL) FREQ UENCIES 
(MODE) SHAPES 


((FOR) STRUCT URE < I abe I > (,) 


'(FOR) LOAD ING < I abe I > ( TIME ) STEP S integer I !st>" 
(FOR) MODE S < integer list> 


Modal loads output request: 

OUTPUT MODAL LOAD S ((FOR) STRUCT URE <label>) (,) 

"(FOR) LOADING < I abel > 

( FOR ) MODES <! nteger I !st> 


q ± User - supplied .Mode . Shapes . 


Command sequence: 


ALTERNATE (MODES) < I abe I > (( TITLE ) <strlng>) 
<specif 1 cation of DOF order: U V W UX ...> 
MOD E <mode number: I nteger> 

C<node number : I nteger > [<DOF val ue: number >33 


Spec! f I cat! on of Dynamic 


Element declaration for intermediate level structure: 

ELEMENT 1 TYPE structure name> CONDENSED (,) 

RETAIN ED ( MODES < i nteger I i st> ) ) 

I NODES <i nteger > ( ( 

USE ALTERNATE ( MOD ES) <label> j 
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Spec! f I cat! on of Initial Conditions 
Command sequence: 

INIT IAL COND ITIONS < i abe I > (( TITLE ) <strlng>) 

s _ 

DISPL ACEMENTS 

C<node l!st> <D0F l!st> - <number>3 




VELOCITIES 

[<node llst> <D0F llst> - <number>3 
USE ( NONLIN EAR) DISPL ACEMENTS ((FOR) STRUCT URE < I abei >) (,) 
( FOR ) LOADING < I abe I > ( STEP <lnteger» 




Speclil cation q± .ilyjmix: Isadiaa CoiiditLaa 

Input of support accelerations as F(x,y,z)s 

LOAD ING < I abe I > (( TITLE ) <strfng>) 

(NODAL) ACCEL ERATIONS 
[<node 1 1 s+> <DOF llst> <number>] 

Definition of loading conditions 

LOADING <label> (( TITLE ) <string>) 

" DYNAM IC 

NONLIN EAR 

Definition of G(t) : 

For a harmonic variation of G(t)s 

HARMONIC PERIOD <number> ( PHASE ( ANGLE ) <number» ( # ) 
( COMBINE ) [ < label > (FACTOR) <number> U ] 

For a general variation of G(+)s 


GENERAL (COMBINE) 

< label > 

"(TIMES ) <n umber lfst> 




\ FACTORS f J 
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For an Impulsive variation of G(+): 

IMPULSIVE (SHAPE) ( HALF_S IN E ) DURAT ION <number> (,) 

) RECTANG ULAR [ 

) P0S_TR I ANGULAR ( 

I NEG_TR I ANGULAR J 

( COMBINE ) C < I abe I > ( FACT OR) <number> 2 


Step size definitions 

C ( TIME ) STEP S <lnteger I Ist> ( ( TITLE ) <strlng>) <number llst> (,) 
(SECONDS) ] 

Definition of results saved In the data base: 

( SAVE ( TIME ) STEP S < Integer llst» 

Note that the last step computed Is always saved, even if not In the 
Integer list or If the command Is not given. 


Spec meat , i on of Transient Analysis 

Definition of the transient analysis methods 

TRANS IENT ANALYSIS (TYPE) ( MODE_SUPER POSITION 

\ NEWMARK 

j CENTRALJ3I FF ERENCE 

I • 

I e 

I ® 

PROPERT I ES <define properties unique to each type> 



Computation requests 


COMPUTE 


NONLIN EAR" 

DYNAMIC 


DISPL ACEMENTS ((FOR) STRUCT URE <Jabel>) 


LOADING < I abe I > ( TIME ) STEP S < Integer I Ist> 
INIT IAL COND ITIONS <label> 

INCLUDE MODES <lnteger I Ist> 




267 


Output request: 


OUTPUT 


" DYNAM IC 

NONLINEAR 


DISPL ACEMENTS 

VELOCITIES 

STRAIN S 

STRESS ES 


(< Integer I ist>) 


(,) 


((FOR) STRUCT URE <label>) (»> 

(FOR) LOAD ING <label> ( TIME ) STEP S < i nteger Hst> 


icaaLLLsatLan fit Shock Sp e ct r um Analysis 


Definition of the spectrum: 

( SHOCK ) SPECTRUM < 1 abe I > (( TITLE ) <strlng>) 

~( DISPL ACEMENTS ) 

< VELOCITIES ? <number I ist> 

( ACCEL ERATIONS ) 

PERIOD S \ <number 1 1 st> 
FREQUENCIES } 

DIRECT IONS (,) 

<node I !st> 

Output requests 


! _U. ) <dl rection coslnesnumber> 

V ' 
w j 


OUTPUT DYNAMIC 


' DISPLA CEMENTS 
VELOCITIES 
STRESS ES 
STRAIN S 

PART I C I PAT I ON_FACTORS 


((FOR) STRUCT URE <label» (») 
(FOR) (SHOCK) SPECTURM <label> 
(FOR) MODES 


(<! nteger I lst>) ( y ) 


< I nteger llst> 
SRSS 

PEAK_SRSS 
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Commands 


Trace command: 


trace ; r nonlin ear - ! ( solut i on ) 

DYNAM IC 

Destroy command: 

DESTROY f NONLIN EAR"] RESULT S ( FOR ) STRUCT URE <!abel> (,) 
DYNAM IC 

( (FOR) LOADING < I abe I > ( TIME ) STEP S < Integer I Ist> 
(FOR) ( SHOCK ) SPECTRUM <label> 

[( FOR ) MOD ES < Integer I tst> 



ooooooooooooooooo 
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7.4 fLa m& l e 

The following sections contain example Input data illustrating the 
use of the foregoing commands for dynamic analysis of some simple plane 
structures. Each example problem Is liberally commented to explain the 
analysis process. The substructured nonlinear analysis of Section 7.4.5 
deserves additional discussion and Is described In detail in that sec- 
tion. 


7.4.1 Standard Linear 


- IXtatlsn An aly sis 


*RUN FINITE 


EXAMPLE INPUT NO. 1 FOR DYNAMIC ANALYSIS 



(REFERENCE FIGURE 7.1) 

THIS EXAMPLE ILLUSTRATES THE INPUT NECESSARY TO DESCRIBE 
A LINEAR, THREE ELEMENT. PLANE FRAME AND TO PERFORM 
A FREQUENCY ANALYSIS OF THE STRUCTURE. IN ANTICIPATION 
OF A TRANSIENT ANALYSIS BY MODE SUPERPOSITION. THE 
LOADS ARE OUTPUT IN MODAL COORDINATES. THE PROBLEM IS 
RESTARTED AND THE TRANSIENT ANALYSIS IS INVOKED WITH ONLY 
SELECTED MODES INCLUDED. INITIAL CONDITIONS ARE ALSO 
DEFINED AND INCORPORATED INTO THE TRANSIENT ANALYSIS. 


STRUCTURE FRAME 
NUMBER OF ELEMENTS 3 NODES 4 

ELEMENTS 1,3 TYPE PLANEFRAME LUMPED E 30000 G 12000 AX 20. , 
AY 5.877 IZ 724 DENSITY 0.00074 
ELEMENT 2 TYPE PLANEFRAME CONSISTENT E 30000 G 12000 AX 14.4, 
AY 3.4 IZ 273. DENSITY 0.00074 

COORDINATES 


1 0.0 

0.0 

2 0.0 

96.0 

3 96.0 

96.0 

4 96.0 

0.0 

INCIDENCES 


1 1 2 




OOO O O OOO OOO OOO OOO 


270 


2 2 3 

3 3 4 

CONSTRAINTS 
1 ,4 ALL = 0.0 

LOAD I NG MOTOR 

ELEMENT LOADS FOR TYPE PLANEFRAME 
2 CONCENTRATED FORCE Y P 300.0 L 48.0 

MASS 

NODAL MASS 
2.3 U V W 1.0 
ELEMENT MASS FOR TYPE PLANEFRAME 

2 UNIFORM U V W W 0.05 $ THE SECOND W INDICATES AN INTENSITY 

USE LOADING MOTOR G 386.4 

STATIC LOADING PATTERN: F(X,Y,Z) 

LOADING PATTERN 

ELEMENT LOADS FOR TYPE PLANEFRAME 
1 LINEAR FORCE Y LA 0.0 LB 1.0 WA 0.0 WB 1.0 

DYNAMIC LOADING CONDITION: G(T) 

LOADING SHAKE 
DYNAMIC 

HARMONIC PERIOD 0.04 PHASE 0.0 COMBINE PATTERN 6.0 
TIME STEPS 1-100 0.0 TO 10.0 BY 0.1 

SAVE STEPS 10-100 BY 10 

SPECIFY FREQUENCY ANALYSIS PARAMETERS 

FREQUENCY ANALYSIS TYPE SUBSPACE 
PROPERTIES 

'CONVERGENCE TOLERANCE 1 .0E-08, 

MAXIMUM ITERATIONS 13, 

SHIFT EVERY 2 MODES, 

MAXIMUM MODES 4, 

RANGE MIN 0.0 MAX 50.0 

REQUEST COMPUTATION OF FREQUENCIES AND MODES 

COMPUTE NATURAL FREQUENCIES MODE SHAPES STRUCTURE FRAME 

OUTPUT FREQUENCIES SHAPES STRUCTURE FRAME MODES 1-4 
OUTPUT MODAL LOADS STRUCTURE FRAME LOADING SHAKE /MODES 1-4 

STOP 


*RUN FINITE FILES = 20,, 22,23 
C 
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RESTART OF THE LINEAR PLANEFRAME TO PERFORM THE 
TRANSIENT ANALYSIS AFTER EXAMINING THE MODAL LOADS. 


ACCESS STRUCTURE FRAME 

SPECIFY TRANSIENT ANALYSIS PARAMETERS 
(USE DEFAULT PROPERTIES) 

TRANSIENT ANALYSIS TYPE MODE_SUPERPOSITION 


DEFINE INITIAL CONDITIONS FOR DISPLACEMENTS. 

(VELOCITIES NOT REQUIRED SINCE WE HAVE NO DAMPING.) 

INITIAL CONDITIONS DEAD_SHAPE 
DISPLACEMENTS 
2,3 U * 0.1 

2 THETAZ » -0.085 

3 THETAZ * 0.085 

COMPUTE DYNAMIC DISPLACEMENTS STRUCTURE FRAME LOADING SHAKE, 

TIME STEPS 1-25 INITIAL CONDITIONS DEAD_SHAPE INCLUDE, 

MODES 1, 3, 4 

OUTPUT DYNAMIC DISPLACEMENTS 2,3 STRUCTURE FRAME LOADING SHAKE, 

STEPS 10, 20, 25 $ RECALL THAT STEP 25 IS SAVED EVEN THOUGH 

I DIDN’T REQUEST IT. 

OUTPUT WIDE BY ELEMENT DYNAMIC STRESSES STRAINS ALL STRUCTURE FRAME, 
LOADING SHAKE STEPS 10, 20, 25 

STOP 
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3 - Node Number Y 

(?)- Element Number f 

NODAL DOF: C • ~ 
^9z 


Figure 7.1 — Three Element, Plane Frame (Example Input #1-3) 
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7.4.2 Standard Linear Structure - .ShO.c. k Spectrum Ana I y . s j s 
*RUN FINITE 


EXAMPLE INPUT NO. 2 FOR DYNAMIC ANALYSIS 


(REFERENCE FIGURE 7.1) 

THIS EXAMPLE ILLUSTRATES THE INPUT NECESSARY TO DESCRIBE 
A LINEAR, THREE ELEMENT, PLANE FRAME AND TO PERFORM 
A SHOCK SPECTRUM ANALYSIS OF THE STRUCTURE. THE SHOCK 
SPECTRUM CAN CONTAIN BOTH HORIZONTAL AND VERTICAL 
COMPONENTS. WE CAN GET THE PARTICIPATION FACTORS 
(ALSO KNOWN AS EFFECTIVE MODAL MASS) PRIOR TO REQUESTING 
THAT SPECTRAL DISPLACEMENTS, STRESSES, AND STRAINS BE 
COMPUTED. TWO METHODS ARE AVAILABLE FOR SUMMING THE 
SPECTRAL VALUES FOR EACH MODE ; SRSS AND PEAK_SR$$. 

SRSS IS THE SQU ARE-ROOT-OF- THE- SUM-OF- THE- SQUARES METHOD. 
PEAK_SRSS TAKES THE MODE WITH THE LARGEST PARTICIPATION 
FACTOR AND ADDS TO THAT THE SRSS OF THE REMAINING MODES 
THAT ARE INCLUDED IN THE ANALYSIS. 


STRUCTURE FRAME 
NUMBER OF ELEMENTS 3 NODES 4 

ELEMENTS 1,3 TYPE PLANEFRAME LUMPED E 30000 G 12000 AX 20. , 

AY 5.877 IZ 724 DENSITY 0.00074 
ELEMENT 2 TYPE PLANEFRAME CONSISTENT E 30000 G 12000 AX 14.4 , 
AY 3.4 IZ 273. DENSITY 0.00074 
C 

COORDINATES 

1 0.0 0.0 

2 0.0 96.0 

3 96.0 96.0 

4 96.0 0.0 

C 

INCIDENCES 
1 1 2 

2 2 3 

3 3 4 

C 

CONSTRAINTS 
1,4 U V * 0.0 

C 

LOADING MOTOR 

ELEMENT LOADS FOR TYPE PLANEFRAME 
2 CONCENTRATED FORCE Y P 300.0 L 48.0 
C 

MASS 



OO 0000000000*000 o o ooo ooooo 


274 


NODAL MASS 
2,3 U V W 1.0 

ELEMENT MASS FOR TYPE PLANEFRAME 
2 UNIFORM U V W W 0.05 
USE LOADING MOTOR G 386.4 

SHOCK SPECTRUM FA I L_SAFE ”5 PERCENT DAMPING" 
VELOCITIES 2 30 30 12 

PERIODS .05 .60 4.5 10. 

DIRECTIONS 1,4 U 0.866 V 0.5 


SPECIFY FREQUENCY ANALYSIS PARAMETERS 

FREQUENCY ANALYSIS TYPE SUBSPACE 
PROPERTIES 

CONVERGENCE TOLERANCE 1 .0E-08, 

MAXIMUM ITERATIONS 13, 

SHIFT EVERY 2 MODES, 

MAXIMUM MODES 4, 

RANGE MIN 0.0 MAX 50.0 

REQUEST COMPUTATION OF FREQUENCIES AND MODES 

COMPUTE NATURAL FREQUENCIES MODE SHAPES STRUCTURE FRAME 

OUTPUT FREQUENCIES SHAPES STRUCTURE FRAME MODES 1-4 
OUTPUT DYNAMIC PARTICIPATION FACTORS STRUCTURE FRAME SHOCK, 
SPECTRUM FA I L_SAFE MODES ALL 


STOP 


RUN FIN 


ITE FILES = 20,, 22 ,23 


RESTART OF THE LINEAR PLANEFRAME TO PERFORM THE 
SHOCK SPECTRUM ANALYSIS AND COMPUTE SPECTRAL STRESS 
AND STRAINS. AT THIS POINT, WE HAVE HAD AN OPPORTUNITY 
TO EXAMINE THE PARTICIPATION FACTORS AND SEE THAT ONLY 
THREE OF THE MODES HAVE ANY SIGNIFICANT CONTRIBUTION TO 
THE SPECTRAL RESPONSE. 


ACCESS STRUCTURE FRAME 


OUTPUT DYNAMIC DISPLACEMENTS STRESSES STRAINS STRUCTURE FRAME, 
SPECTRUM FAIL SAFE MODES 1-3 SRSS PEAK SRSS 

STOP 
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7.4.3 S tandard Nonlinear Structure 


*RUN FINITE 


EXAMPLE INPUT NO. 3 FOR DYNAMIC ANALYSIS 


(REFERENCE FIGURE 7.1) 

THIS EXAMPLE ILLUSTRATES THE INPUT NECESSARY TO DESCRIBE 
A NONLINEAR, THREE ELEMENT, PLANE FRAME AND TO PERFORM 
A TRANSIENT ANALYSIS OF THE STRUCTURE. AT SOME SELECTED 
TIME STEP, THE TRANSIENT ANALYSIS IS SUSPENDED AND A 
FREQUENCY ANALYSIS IS PERFORMED WITH THE CURRENT STRUCTURE 
STIFFNESS AND MASS. THE GENERALIZED NEWMARK OPERATOR IS 
USED FOR THE TIME-HTORY INTEGRATION. THIS OPERATOR 
CAN BE USED AS EITHER AN EXPLICIT OR IMPLICIT INTEGRATOR 
AND HAS THE ABILITY TO CONTROL SPURIOUS DAMPING. 


MATERIAL STEEL TYPE VON._M!SES 
PROPERTIES S IGNAL_YIELD 
USE STRESS- STRAIN FUNCTION SEGMENTAL 

PROPERTIES E 30000. NU 0.3 STRA I N__HARDEN I NG TENSYIELD 30, 
COMPYIELD 30. TENS 1 0N„SL0PE 3000. COMP„$LOPE 3000. 


STRUCTURE FRAME 
NUMBER OF ELEMENTS 3 NODES 4 

ELEMENTS 1,3 TYPE PLANEFRAME LUMPED E 30000 G 12000 AX 20. , 
AY 5.877 IZ 724 DENSITY 0.00074 
ELEMENT 2 TYPE PLANEFRAME CONSISTENT MATERIAL STEEL E 30000 , 
G 12000 AX 14.4 AY 3.4 IZ 273. DENSITY 0.00074 
C 

COORDINATES 

1 0.0 0.0 

2 0.0 96.0 

3 96.0 96.0 

4 96.0 0.0 

C 

INCIDENCES 
1 1 2 

2 2 3 

3 3 4 

C 

CONSTRAINTS 
1 .4 ALL = 0.0 
C 

LOADING MOTOR 

ELEMENT LOADS FOR TYPE PLANEFRAME 
2 CONCENTRATED FORCE Y P 300.0 L 48.0 
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MASS 

NODAL MASS 
2,3 U V W 1.0 
ELEMENT MASS FOR TYPE PLANEFRAME 
2 UNIFORM U V W W 0.05 
USE LOADING MOTOR G 386.4 

LOADING PATTERN 

ELEMENT LOADS FOR TYPE PLANEFRAME 
1 LINEAR FORCE Y LA 0.0 LB 1.0 WA 0.0 WB 1.0 

LOADING SHAKE 
NONLINEAR DYNAMIC 

HARMONIC PERIOD 0.04 PHASE 0.0 COMBINE PATTERN 6.0 
TIME STEPS 1-100 0.0 TO 10.0 BY 0.1 

SAVE STEPS 10-100 BY 10 

SPECIFICATION OF NONLINEAR SOLUTION PARAMETERS 

TRACE NONLINEAR SOLUTION 

CONVERGENCE TEST NORM RESIDUAL LOADS TOLER 1.5 INCLUDE TOTAL, 
REACTIONS 

UPDATE STIFFNESS EVERY STEP 
TERMINATE IF NONCONVERGENT 

SPECIFICATION OF TRANSIENT ANALYSIS SOLUTION PARAMETERS 

TRANSIENT ANALYSIS TYPE NEWMARK 

PROPERTIES ALPHA 0.0 BETA 0.5 GAMMA 0.25 

REQUESTS FOR COMPUTATION AND OUTPUT. 

COMPUTE NONLINEAR DYNAMIC DISPLACEMENTS STRUCTURE FRAME LOADING, 
SHAKE TIME STEPS 1-30 

OUTPUT DYNAMIC NONLINEAR STRESSES 2 STRUCTURE FRAME LOADING SHAKE, 
STEPS 10.20.30 


STOP 


UJN FINITE FILES » 20,, 22 ,23 


RESTART OF NONLINEAR PLANE FRAME TO PERFORM THE 
FREQUENCY ANALYSIS AND GET THE NATURAL FREQUENCIES 
AND MODE SHAPES. 


ACCESS STRUCTURE FRAME 


SPECIFICATION OF FREQUENCY ANALYSIS 
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FREQUENCY ANALYSIS TYPE JACOBI 
PROPERTIES 
MAX_SWEEPS 15, 

CONVERGENCE TOLERANCE 1 .0E-08 


REQUEST COMPUTATION AND OUTPUT OF FREQUENCIES AND MODES 

COMPUTE NONLINEAR NATURAL FREQUENCIES MODE SHAPES STRUCTURE FRAME, 
LOADING SHAKE TIME STEP 20 

OUTPUT NONLINEAR FREQUENCIES SHAPES STRUCTURE FRAME MODES ALL, 
LOADING SHAKE STEP 20 


STOP 


AT THIS POINT, WE CAN "ACCESS STRUCTURE..." AGAIN 
AND CONTINUE WITH THE TRANSIENT ANALYSIS 
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7.4.4 Substructured Linear Analysis 


*RIJN FINITE 


EXAMPLE INPUT NO. 4 FOR DYNAMIC ANALYSIS 


(REFERENCE FIGURE 7.2) 

THIS EXAMPLE ILLUSTRATES THE INPUT NECESSARY TO DESCRIBE 
A LINEAR, MULTILEVEL SUBSTRUCTURED MODEL AND TO PERFORM 
A FREQUENCY ANALYSIS OF THE HIGHEST LEVEL STRUCTURE. 
FEATURES OF DYNAMICS ILLUSTRATED ARE GUYAN REDUCTION AND 
MODAL SYNTHESIS FOR CONDENSING THE SUBSTRUCTURES AND THE 
METHOD FOR CARRYING FORWARD SUBSTRUCTURE MASS (IE. NO 
SPECIAL CONSIDERATION IS GIVEN TO BRINGING UP MASS). 

THE ENTIRE STRUCTURAL SYSTEM IS BUILT OUT OF ONE STAND- 
ALONE ELEMENT. A PLANE FRAME ELEMENT. 


ELEMENT BAR TYPE PLANEFRAME CONSISTENT E 30000. G 12000. AX 20.0, 

AY 5.877 IZ 724 DENSITY 0.00074 

COORDINATES 

1 0.0 0.0 

2 96.0 0.0 

DEFINE LOWEST LEVEL STRUCTURE, A TRUSS WITH THREE BAYS. 

STRUCTURE THREE_BAY 
NUMBER OF ELEMENTS 7 NODES 5 
ELEMENTS ALL TYPE BAR 

1, 6, 7 ROTATION SUPPRESSED 

2, 4 ROTATION Z 60. 

3, 5 ROTATION Z -60. 

INCIDENCES 
1 1 2 

2 3 1 

3 1 4 

4 4 2 

5 2 5 

6 3 4 

7 4 5 

ADD MASS TO ONE OF THE LOWER CHORDS. 

MASS 

ELEMENT MASS FOR TYPE PLANEFRAME 
7 UNIFORM U V W W 0.0003 
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CONDENSE OUT NODE 4 USING GUYAN REDUCTION 

STRUCTURE THREE_C0N 
NUMBER OF ELEMENTS 1 NODES 4 
ELEMENT 1 TYPE THREE_BAY CONDENSED 

INCIDENCES 
1 3 15 2 


STICK TWO OF THE FRAMES TOGETHER AND CLOSE THE GAP 
AT THE TOP WITH A BAR ELEMENT. 

STRUCTURE SPAN 

NUMBER OF ELEMENTS 3 NODES 7 
ELEMENTS 

1, 3 TYPE THREE_C0N ROTATION SUPPRESSED 
2 TYPE BAR ROTATION SUPPRESSED 

INCIDENCES 

1 12 3 4 

2 4 5 

3 3 5 6 7 

ADD A LITTLE MORE MASS TO SOME SELECTED NODES 


MASS 

NODAL MASS 

4, 5 U V W 0.5 

DEFINE FREQUENCY ANALYSIS PARAMETERS FOR THIS STRUCTURE 
SINCE IT WILL BE CONDENSED USING MODAL SYNTHESIS 

FREQUENCY ANALYSIS TYPE SUBSPACE 
PROPERT I ES 

CONVERGENCE TOLER 1.0E-08, 

MAX ITERATIONS 10, 

MAX MODES 5 


CONDENSE SPAN VIA MODAL SYNTHESIS 

STRUCTURE SPAN_CON 

NUMBER OF ELEMENTS 1 NODES 2 

ELEMENT 1 TYPE SPAN CONDENSED RETAINED MODES 5 


INCIDENCES 
1 1 6 


BUILD THE TWO SPAN BRIDGE 

NOTE THAT THE -mETAZ DOF AT THE MIDDLE SUPPORT KEEPS 
THIS STRUCTURE FROM BECOMING TWO SIMPLE SPANS. 
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ALSO NOTE THAT ALTHOUGH THIS STRUCTURE ONLY HAS 3 NODES 
WITH 3 GEOMETRIC DOF EACH, IT HAS 10 GENERALIZED DOF FROM 
THE RETAINED NORMAL MODES FOR A TOTAL OF 19 DOF. 

STRUCTURE TV/0_SPAN 
NUMBER OF ELEMENTS 2 NODES 3 
ELEMENTS ALL TYPE SPAN_CON ROTATION SUPPRESSED 
C 

INCIDENCES 
1 1 2 

2 2 3 

C 

CONSTRAINTS 
1 U V = 0.0 
23 V - 0.0 
C 

FREQUENCY ANALYSIS TYPE SUBSPACE 
PROPERTIES 

CONVERGENCE TOLER 1 .0E-08, 

MAX ITERATIONS 10, 

MAX MODES 10 
C 

COMPUTE FREQUENCIES STRUCTURE TWO_SPAN 
C 

OUTPUT FREQUENCIES SHAPES STRUCTURE TWO_SPAN MODES ALL 
C 


STOP 




7.4.5 Substructured. Non I Inear 


The Input listing In this section describes a simple, two-span 
bridge. Nonlinearity is introduced Into the example by addition of a 
nonlinear bar element over the center support. The example uses six 
levels of substructures with the nonlinear element added at the highest 
level. To facilitate the description of the input sequence, line 
numbers are placed before each FINITE command line- The numbers are not 
part of the commands. They serve only as reference numbers. Comment 
lines are not numbered. The structure is illustrated in Figure 7.3. 

Lines 2-6 define a stand-alone element, which is used to construct 
the majority of the final structure. Lines 7-23 describe the lowest 
level substructure called THREEJ3AY. No frequency analysis parameters 
are defined since this structure will be condensed using Guyan reduc- 
tion. The condensed version of THREEJ3AY Is named structure PIECE. No 
additional input is required to Invoke the condensation process; Guyan 
reduction Is the default procedure adopted by the system. 

Lines 29-48 describe structure HALF which contains two copies of 
substructure PIECE and one copy of stand-alone element BAR. To il- 
lustrate its use In substructures, secondary mass is applied in lines 
38-40. A frequency analysis method Is defined, lines 44-48, so the 
structure can be condensed by modal synthesis. Structure HALF_C0N. 
lines 49-56, is the condensed version of structure HALF. It is neces- 
sary to carry forward the loads from HALF but the mass Is automatically 
Included In the reduction process. 
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A nonlinear material model is defined in lines 57-61. This model 
is required for the nonlinear element used in structure BRIDGE. This 
highest level structure is composed of two condensed substructures, 
HALF_C0N. and one simple element, TYPE PLANEFRAME (see lines 62-67). 
Damping and dynamic loading are defined in lines 78-88 with the tran- 
sient analysis specification and requests for computation and output in 
I ines 89-96. 


1 *RUN FINITE 
C 

C 

C EXAMPLE INPUT NO. 5 FOR DYNAMIC ANALYSIS 

C ======S==S=================== * ===:==SS3:==:= 

C 

C (REFERENCE FIGURE 7.3) 

C 

C THIS EXAMPLE ILLUSTRATES THE INPUT NECESSARY TO DESCRIBE 

C A NONLINEAR, MULTILEVEL SUBSTRUCTURED MODEL AND TO PERFORM 

C A TRANSIENT ANALYSIS OF THE STRUCTURE WHEN IT IS SUBJECTED 

C TO A GENERAL DYNAMIC LOADING. THE NONLINEARITY IS 

C RESTRICTED TO A MATERIALLY NONLINEAR ELEMENT AT THE HIGHEST 

C LEVEL. THE STRUCTURE IS THE TRUSS FROM EXAMPLE NO. 4 WITH 

C THE GAP OVER THE CENTER SUPPORT CLOSED BY THE NONLINEAR 

C ELEMENT. 

C 

C 

2 ELEMENT BAR TYPE PLANEFRAME CONSISTENT E 30000. G 12000. AX 20.0, 

3 AY 5.877 IZ 724 DENSITY 0.00074 
C 

4 COORDINATES 

5 1 0.0 0.0 

6 2 96.0 0.0 

C 

C DEFINE LOWEST LEVEL STRUCTURE, A TRUSS WITH THREE BAYS. 

C 

7 STRUCTURE THREE_BAY 

8 NUMBER OF ELEMENTS 7 NODES 5 

9 ELEMENTS ALL TYPE BAR 


10 

1, 6, 7 

ROTATION 

SUPPRESSED 

11 

2, 4 

ROTATION 

Z 60. 

12 

3, 5 

ROTATION 

Z -60. 


C 



13 

INCIDENCES 



14 

1 1 2 



15 

2 3 1 



16 

3 1 4 
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17 4 4 2 

18 5 2 5 

19 6 3 4 

20 7 4 5 

C 

C ADD MASS TO ONE OF THE LOWER CHORDS. 

C 

21 MASS 

22 ELEMENT MASS FOR TYPE PLANEFRAME 

23 7 UNIFORM U V W W 0.0003 
C 

C 

C CONDENSE OUT NODE 4 USING GUYAN REDUCTION 

C 

24 STRUCTURE PIECE 

25 NUMBER OF ELEMENTS 1 NODES 4 

26 ELEMENT 1 TYPE THREE_BAY CONDENSED 
C 

27 INCIDENCES 

28 1 3152 

C 

C 

C STICK TWO OF THE FRAMES TOGETHER AND CLOSE THE GAP 

C AT THE TOP WITH A BAR ELEMENT. 

C 

29 STRUCTURE HALF 

30 NUMBER OF ELEMENTS 3 NODES 7 

31 ELEMENTS 

32 1,3 TYPE PIECE ROTATION SUPPRESSED 

33 2 TYPE BAR ROTATION SUPPRESSED 

C 

34 INCIDENCES 

35 112 3 4 

36 2 4 5 

37 3 3 5 6 7 

C 

C ADD A LITTLE MORE MASS TO SOME SELECTED NODES 

C 

38 MASS 

39 NODAL MASS 

40 4, 5 U V W 0.5 
C 

C ADD THE PATTERN OF LOAD TO BE USED IN THE DYNAMIC 

C LOADING 

C 

41 LOADING CENTER_SPAN 

42 NODAL LOADS 


43 3 FORCE Y P -1.0 
C 

C 

C DEFINE FREQUENCY ANALYSIS PARAMETERS FOR THIS STRUCTURE 

C SINCE IT WILL BE CONDENSED USING MODAL SYNTHESIS 

C 

44 FREQUENCY ANALYSIS TYPE SUBSPACE 
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45 

46 

47 

48 


49 

50 

51 

52 

53 


54 

55 

56 


57 

58 

59 

60 
61 


62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 


C 

C 

C 

C 


C 


C 

C 

C 


C 

C 

C 

C 

c 


c 

c 

c 

c 

c 

c 

c 


c 


c 


c 

c 


PROPERTIES 

CONVERGENCE TOLER 1.0E-08, 
MAX ITERATIONS 10. 

MAX MODES 5 


CONDENSE HALF VIA MODAL SYNTHESIS 

STRUCTURE HALF_CON 

NUMBER OF ELEMENTS 1 NODES 4 

ELEMENT 1 TYPE HALF CONDENSED RETAINED MODES 3 


INCIDENCES 
1 12 6 7 


CARRY FORWARD THE LOADS FROM HALF 

LOADING CENTER_CON 
EXTERNAL ELEMENT LOADS 
1 CENTS R_SPAN 1 .0 


DEFINE THE NONLINEAR MATERIAL 

MATERIAL STEEL TYPE VONJMISES 
PROPERTIES SIGNAL YIELD 
USE STRESS” STRAIN FUNCTION SEGMENTAL 
PROPERTIES E 1 NU 0 STRAIN HARDENING, 
TENSYIELD 1 TENS I ON_SLOPE .1 


BUILD THE TWO SPAN BRIDGE 

CLOSE THE GAP OVER THE CENTER SUPPORT WITH A 
NONLINEAR BAR. 

STRUCTURE BRIDGE 
NUMBER OF ELEMENTS 3 NODES 7 
ELEMENTS 

1, 3 TYPE HALF_CON ROTATION SUPPRESSED 
2 TYPE PLANEFRAME MATERIAL STEEL CONSISTENT AX 20.0 
AY 5.877 IZ 724 DENSITY 0.00074 

COORDINATES 

4 0.0 0.0 

5 96.0 0.0 

INCIDENCES 

1 12 3 4 

2 4 5 

3 3 5 6 7 


9 
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75 CONSTRAINTS 

76 1 U V = 0.0 

77 36 V = 0.0 
C 

C APPLY DAMPING TO THE STRUCTURE BY USING RAYLEIGH DAMPING 

C 

78 DAMPING RAYLEIGH FREQUENCIES 2.0 12.0 PERCENTS 1.0 3.4 

C 

C 

C APPLY THE LOAD PATTERN TO EACH SPAN 

C 

79 LOADING PATTERN 

80 EXTERNAL ELEMENT LOADS 

81 1,3 CENTER_CON 1.0 

C 

C DEFINE THE DYNAMIC, NONLINEAR LOAD. 

C 

82 LOADING SHAKE 

83 DYNAMIC NONLINEAR 
C 

84 GENERAL COMBINE PATTERN, 

85 FACTORS 0.0 100. 40. -50. 100. 0. 50. -90. 0.0, 

86 TIMES 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 

C 

87 TIME STEPS 1-100 0.0 TO 1.0 BY .01 

88 SAVE STEPS 5-100 BY 5 

C 

C DEFINE THE TRANSIENT ANALYSIS 

C 

89 TRANSIENT ANALYSIS TYPE NEWMARK 

90 PROPERTIES 

91 ALPHA 0.0, 

92 BETA 0.5, 

93 GAMMA 0.25 

C 

94 COMPUTE NONLINEAR DYNAMIC DISPLACEMENTS STRUCTURE BRIDGE, 

95 LOADING SHAKE TIME STEPS 1-25 

C 

96 OUTPUT NONLINEAR DISPLACEMENTS LOADING SHAKE TIME STEPS 5-25 BY 5 
C 

c 

97 STOP 
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CHAPTER 8 


SUMMARY AND TOPICS FOR FURTHER STUDY 

8.1 Summary 

Comprehensive dynamic analysis of complex structural systems by the 
finite element method can be an expensive, if not impossible under- 
taking. Existing software systems capable of achieving some economy 
suffer a limited scope. The need exists for a general purpose FEM 
system which is capable of dynamic analysis of arbitrary structures. 
This capability includes structures experiencing geometric and material 
nonlinearities. In order to achieve an economical solution, multilevel 
substructur ing is seen as a requisite modeling approach. It is the pur- 
pose of this work to bring together the individual, isolated topics of 
multilevel substructured modeling, dynamic analysis by the FeM, and non- 
linear continuum mechanics into the design of a comprehensive, general 
purpose, finite element package. The resulting software will be used to 
perform numerical experiments to explore the behavior of the proposed 
modal synthesis technique in a multilevel substructured environment. 
The factors studied will include the economics, accuracy, and analyst 
interaction required to perform modal synthesis. 

Implementation of multilevel substructur ing for static analysis of 
linear and nonlinear structures has been discussed In detail. The suc- 
cess of the effort Is dependent upon the schemes used for data storage 
and retrieval, equation solving, and user definition of the model. It 
was shown that static results are equivalent for both substructured and 
standard models. Economy in the solution via substructured modeling was 
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demonstrated in two numerical examples. 

Dynamic reduction of the stiffness and mass matrices has been iden- 
tified as the pivotal process in accurately representing complex struc- 
tures as simplified models for dynamic analysis. A number of the 
various methods currently available for dynamic reduction have been 
identified by a review of the open literature. Guyan reduction and the 
f ixed- i nterface method have been chosen for Incorporation into the 
general purpose FEM software system. 

Eigenproblem solution and transient response analysis are the most 
computationally expensive operations In the dynamic analysis of struc- 
tural systems. Their proper Implementation and use Is essential to the 
success of the dynamic analysis. A brief review of these processes and 
an examination of their use in a multilevel substructured environment 
was given. The most effective eigenproblem solution methods have been 
identified while transient response analysis was discussed in more 
general terms. 

Using matrix notation, the nonlinear equations of continuum 
mechanics were derived. Two formulations. Total Lagrangian (T. L.) and 
Updated Lagrangian (U. L.) were described In detail. Both formulations 
were shown to derive from a common definition of the rate of work per 
unit mass and thus should provide identical analysis results. Dif- 
ferences in the computational efficiency of the two formulations were 
shown to arise In the stress rate transformations and in the complexity 
of the nonlinear strain-displacement relations. It was concluded that 
the T. L. formulation has a slight advantage in that no question arises 
regarding the significance of certain nonlinear terms, i.e., all non- 
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linear terms must be included in the formulation. 

The details of a transient solution procedure for a substructured 
nonlinear model based upon an implicit integration operator were 
presented. An implicit scheme was recommended to support dynamic 
analysis since a static solution procedure can be obtained as the 
degenerate case of dynamic analysis. Details of the elemental stiffness 
matrices were derived for both the T. L. and U. L. formulations. 
Specific matrices were listed for the general 2-D case. Qualitative 
comparisons of computational efficiency were made and a T. L» approach 
was recommended for a general software system. The current absence of 
computational evidence in the I Iterature regarding the performance of a 
finite element solution based on each approach does not enable the 
superior approach to be identified. However, an U. L. approach can be 
easily embedded within a T. L. software system. A T. L. approach can- 
not be as easily incorporated into a U. L. based system. Thus, the 
choice of T. L. provides some flexibility for future modifications. 

The POLO-FINITE input language has been extended to include the 
computational features recommended in this report for general purpose 
dynamic analysis. Wherever possible, consistency has been maintained in 
the philosophy and method of defining the substructured model. The com- 
plete command structure was detailed and examples of Its use were 
presented. 
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8.2 Topics for Further Study 

With the definition of the basic requirements for general dynamic 
analysis now available, efforts can be directed to software design. Im- 
plementation, and verification. 

The next task to be performed is the design of a prototype software 
system with the specific goal of demonstrating the applicability of mul- 
tilevel substructur ing in nonlinear dynamic analysis. During the 
literature review, no evidence was found of this having been attempted 
at any level of sophistication. Additional software design topics in- 
clude design of the data structures and processing modules necessary for 
performing the analysis and specification of the formats for convenient 
and selective output of results. 

Later activities include implementation and testing of the system 
In the POLO-FINITE structural mechanics software system. The perfor- 
mance of the system will be evaluated over a broad range of structural 
types including general substructure geometry and I I near/non I i near 


response. 
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